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Methods  for  enhancing  the-  parameter  estimation  process 

/ 

of  phptrmacokinetic  models  were  examined  and  implemented  for 
-AP^RC;  The  Open  Chamber  Hexane  Model  was  the  primary  focus. 
Optimization  alg  >rithms  were  developed  and  used  to  estimate 
the  unknown  parameters  in  order  that  the  model  match  closely 
with  the  experimental  data.  A  sensitivity  analysis  was  then 
conducted  in  two  parts.  First,  the  sensitivity  of  the  estimated 
parameters  to  the  error  criterion  was  examined.  Second,  the 
use  of  sensitivity  functions  in  designing  future  experiments 
was  demonstrated.  Finally,  recommendations  were  provided 
for  AFAMRL  concerning  future  attempts  in  parameter  estimation. 


I.  INTRODUCTION 


At  the  Air  Force  Aerospace  Medical  Research  Laboratory 
(AFAMRL),  researchers  are  studying  the  effects  of  toxic  substances 
on  mammals.  One  reason  is  that  Air  Force  personnel  are  occasionally 
exposed  to  toxic  substances.  By  understanding  how  these  substances 
affect  man,  the  researchers  can  aid  Air  Force  decision  makers 
in  choosing  those  materials  safest  for  Air  Force  use. 

In  their  work,  the  researchers  generally  follow  three  steps. 
First,  they  develop  a  mathematical  model  that  represents  the 
chemical  processes  in  the  mammal.  This  model  is  referred  to 
as  a  pharmacokinetic  model.  Next,  the  researchers  design 
experiments  to  collect  data.  Finally^  they  estimate  any  unknown 
parameters  in  the  model  by  selecting  the  ones  that  provide  the 
closest  agreement  between  the  model  and  the  data.  With  the 
estimates  of  the  parameters,  the  model  is  completely  specified. 

This  model  is  then  used  to  predict  the  concentration  level  of 
toxic  substances  throughout  the  body. 

The  process  of  developing  the  pharmacokinetic  model  and 
designing  the  experiment  has  been  successfully  done  at  several 


places,  including  W right  Patterson  AFB.  However,  the  process  of 
choosing  the  estimates  requires  limited  resources  to  solve  the  model 
and  compare  these  solutions  with  the  experimental  data.  Thus,  parameter 
estimation  is  generally  confined  to  facilities  that  can  provide 
adequate  support  for  solving  pharmacokinetic  models. 

Before  May  of  1982,  AFAMEL  scientists  conducted  at  of 
their  pharmacokinetic  calculations  at  the  facilitie  f  Dow 

Chemical  Company.  The  facilities  were  beneficial  to  „  research, 

but  it  failed  to  meet  some  of  the  scientists'  needs.  In  addition, 
the  facilities'  location  away  from  Wright  Patterson  AFB  proved 
to  be  an  inconvenience  to  the  scientists.  For  these  reasons, 
two  scientists.  Dr.  Anderson  and  Major  Clewell  looked  to  the 
Mathematics  Department  at  AFIT  to  help  develop  local  support  for 
their  work. 

In  May  of  1982,  Mayberry  (GCS-82D)  initiated  a  thesis  effort 
that  successfully  provided  solutions  of  some  pharmacokinetic  models. 
This  effort  was  the  first  step  towards  in-house  estimation  of 
the  model  parameters  at  Wright  Patterson  AFB.  With  Mayberry's 
computer  programs,  the  researchers  could  try  values  for  the 
unknown  parameters  and  plot  the  model's  results.  By  comparing 
these  results  with  the  actual  data,  the  researchers  would  then 
adjust  the  values  to  make  the  model  agree  more  closely  with  the 
data.  This  trial  and  error  process  of  adjusting  the  values  continued 


until  the  model  matched  closely  with  the  data 


Unfortunately,  this  process  of  parameter  estimation  is  very 
tedious  and  time  consuming.  In  the  past,  researchers  have  typically 
taken  about  one  month  to  estimate  the  parameters.  Yet,  from  their 
work  they  often  discover  changes  necessary  to  the  model.  Then 
they  once  again  must  estimate  the  new  model's  parameters.  In 
essence,  the  parameter  estimation  process  is  almost  never  ending. 
Because  of  this  tedium,  the  researchers  requested  that  a  computer 
be  used  to  free  up  their  time  and  provide  results  much  more  quickly. 

In  addition  to  the  tedium  of  parameter  estimation,  another 
problem  encouraged  the  researchers  to  use  the  computer  for 
parameter  estimation.  Some  pharmacokinetic  models  have  grown 
in  complexity,  in  the  number  of  unknown  parameters,  ana  in  the 
number  of  measurements.  As  a  result,  it  has  become  nearly  impossible 
for  the  researchers  to  vary  several  parameters  manually  in  order  to 
fit  the  model  to  several  data  curves.  For  these  models,  computer 
support  is  almost  necessary  to  estimate  the  parameters. 

Problem  Statement 

Methods  must  be  developed  that  improve  the  parameter  estimation 
process  of  a  pharmacokinetic  model.  In  general,  this  includes  any 
effort  that  helps  speed  the  estimation  process  and  guides  the 
design  of  experiments.  However,  AFAMRL  requires  that  these  methods 


be  used  to  estimate  the  parameters  of  a  model  currently  of  interest 


Scope 

This  thesis  effort  only  deals  with  the  "Open  Chamber  Hexane 
Model".  However,  the  principles  involved  will  be  general  enough 
to  apply  to  most  of  the  models  in  which  AFAMRL  is  interested. 
The  computer  programs  developed  are  mainly  intended  to  estimate 
the  parameters  of  the  Open  Chamber  Hexane  Model.  However,  with 
minor  modifications,  they  may  be  used  to  solve  a  large  number  of 
related  models. 

Assumptions 

It  is  important  to  note  that  the  model  developed  by  AFAMRL 
is  new.  As  a  result,  this  thesis  effort  has  run  parallel  with 
AFAMRL* a  continual  updating  of  the  model.  With  an  eye  focused 
towards  supporting  their  research,  the  author  has  kept  close 
contact  with  Dr.  Anderson  and  Major  Clewell.  Any  changes  made 
by  them  have  been  incorporated  in  this  thesis  effort  when 
possible.  It  is  assumed,  however,  that  the  model  reported  in 
this  thesis  is  the  most  correct  form  known. 

General  Approach 

The  first  step  in  estimating  the  parameters  of  the  Open 


Chamber  Hexane  Model  was  to  get  the  model  and  the  experimental 
data*  With  the  model  and  data,  much  of  the  literature  was  reviewed 
in  order  to  locate  topics  that  might  aid  in  solving  the  problem. 

With  a  weighted  least  squares  criterion  for  parameter  estimation, 
optimization  routines  were  developed  for  estimating  the  parameters. 
The  optimization  routines  were  then  extensively  used  in  a  process 
of  seeking  the  global  minimum*  Finally,  a  sensitivity  analysis 
was  conducted  to  evaluate  the  estimates  and  to  provide  clues 
about  designing  future  experiments. 

Sequence  of  Presentation 

Chapter  H  first  introduces  the  compartments  associated  with  the 
Hexane  Model.  The  differential  equations  are  then  derived 
for  the  model.  The  unknown  parameters  are  specified  next. 

Finally,  the  experimental  data  is  obtained  for  the  purpose 
of  estimating  the  unknown  parameters. 

Chapter  III  reviews  the  literature  in  an  attempt  to  gather 
together  the  background  material  relevant  to  parameter  estimation. 
This  chapter  is  meant  to  be  a  resource  for  this  thesis  effort 
as  well  as  any  future  parameter  estimation  efforts. 

Chapter  IV  derives  the  equations  used  for  optimization  with 


a  weighted  least  squares  error  function.  The  Steepest  Descent  and 
the  Scalar  Gauss-Newton  methods  are  the  algorithms  selected  for 


further  use 


Chapter  V  uses  the  optimization  routines  to  estimate  the 
parameters.  The  final  estimates  and  corresponding  data  plots 
were  then  reviewed  by  AFAMRL.  The  estimates  were  then  accepted 
by  AFAMRL  for  their  applications. 

Chapter  VI  looks  at  what  happens  to  the  error  function 
when  the  final  estimates  are  varied.  This  sensitivity  analysis 
gives  indications  about  the  influence  that  the  estimates  have 
on  the  error  function.  Then,  the  sensitivity  analysis  is 
extended  to  cover  some  aspects  that  are  relevant  in  the  design 
of  future  experiments. 


Ha  Formulation  and  Experimentation:  The  Open  Chamber  Hexane  Model 

This  chapter  specifies  the  items  necessary  before  any  attempt 
can  be  mate  to  estimate  parameters  of  a  pharmacokinetic  model. 

The  first  item  is  the  equations  of  the  pharmacokinetic  model  of 
interest.  The  model  used  in  this  thesis  effort  is  the 
Open  Chamber  Hexane  Model.  The  next  items  to  be  considered  are 
the  unknown  parameters  of  the  model.  Finally,  the  experimental 
data  used  to  estimate  the  unknown  parameters  is  presented. 

The  Open  Chamber  Hexane  Model 

The  chemical  hexane  is  a  substance  contained  in  many  materials 
fuels  widely  used  throughout  the  Air  Force.  As  a  result, 
the  Air  Force  is  interested  in  its  effects  on  man.  The  Open 
Chamber  Hexane  Model  was  specifically  formulated  to  describe 
these  effects.  Currently,  the  Hexane  Model  used  by  AFAMRL 
represents  laboratory  mammals  rather  than  man.  However,  a  few 
changes  to  the  model's  parameters  are  all  that  is  needed  for  it 
to  represent  man.  One  probles  with  the  motel  is  that  not  all 
of  its  parameters  are  known.  Intimating  these  parameters  is 


important  in  order  to  use  the  Hexane  Model  for  Air  Force 
applications. 

To  estimate  some  of  these  parameters,  researchers  first  place 
laboratory  mammals  in  test  chambers.  They  then  fill  the  chamber 
with  hexane  gas.  At  various  instances  in  time,  the  gas  concentration 
level  is  recorded.  Simultaneously,  researchers  record  the 
concentration  levels  in  various  parts  of  the  mammal,  such  as 
the  blood  and  body  water.  Together,  these  measurements  provide 
the  data  set  needed  for  parameter  estimation. 

For  the  experiment,  researchers  may  choose  to  either  control 
the  gas  concentration  level  throughout  time  or  just  set  the 
initial  gas  concentration  level.  In  the  first  case,  the  "Open 
Chamber"  Model  is  used.  The  chamber  is  open  in  the  sense  that 
researchers  maintain  the  desired  input  of  hexane  gas.  The 
second  case  is  known  as  the  "Closed  Chamber"  Model.  Once  the 
experiment  starts,  no  exchange  of  air  is  allowed  between  the 
chamber  and  its  outside  environment.  Thus,  the  gas  concentration 
level  in  the  Closed  Chamber  Model  is  controlled  only  by  the  mammals’ 
own  metabolic  processes  and  not  by  the  experimenter. 


formulate  pharmacokinetic  models,  researchers  first  represent 
the  "“»"■« as  a  set  of  compartments.  Each  compartment  is  typically 
an  organ,  tissue,  or  any  part  of  the  body  where  "the  drug 


concentration  is  assumed  to  be  uniform.”  (MAYB82b:  2.6)  Also, 
separate  compartments  are  used  to  model  different  drugs  that  are 
found  in  the  same  area.  Next,  researchers  connect  these  compartments 
to  show  where  the  chemicals  are  exchanged.  A  more  detailed 
discussion  on  compartmental  modelling  is  provided  by  Mayberry  (MAY82b) 
In  this  section,  the  Open  Chamber  Hexane  Model  will  be  described 
with  compartments. 

When  inhaled,  the  hexane  gas  basically  undergoes  two  types 
of  processes:  (1)  the  distribution  process  and  (2)  the  metabolic 
process.  A  distribution  process  is  the  transfer  of  the  substance 
from  one  location  in  the  body  to  another.  The  transfer  does 
not  change  the  chemical  properties  of  the  substance.  For  example, 
the  transfer  of  the  hexane  gas  from  the  lung  to  the  blood  is  a 
distribution  process.  A  metabolic  process,  on  the  other  hand, 
is  a  process  that  changes  the  substance's  chemical  properties. 

New  chemicals  often  result  as  by-products  during  the  metabolism. 

These  new  substances  are  called  metabolites.  An  example  of 
metabolism  is  the  oxidation  of  a  substance  by  the  liver.  For 
the  Hexane  Model,  the  distribution  process  will  be  discussed 
first. 

In  the  Hexane  Model,  five  compartments  are  primarily  involved 
in  the  distribution  process.  These  compartments  shown  in  Figure  2.1 
are  the  Hexane  in  Lung,  Hexane  in  Rapidly  Perfused  Tissue, 


Hexane  in  Slowly  Perfuse i  Tissue*  Hexane  in  Fat,  and  Hexane  in 
Liver.  The  distribution  process  begins  when  the  inhales 

the  hexane  gas  into  the  lungs.  Froa  the  lungs,  hexane  is  transferred 


to  the  blood.  The  blood  nay  be  considered  as  a  compartment, 
although  this  is  not  shown  in  Figure  2.1.  The  reason  is  that  the 
gas  enters  and  leaves  the  blood  at  a  much  higher  rate  than  the 
other  compartments.  Thus,  the  concentration  level  in  the  blood 
is  just  a  weighted  combination  of  the  concentration  levels  of 
each  of  the  distribution  compartments  in  its  path.  In  essence, 
the  blood  performs  the  function  of  a  transport  medium.  It 
distributes  the  hexane  from  the  lung  compartment  to  the  other 
four  compartments  in  the  model. 

The  other  four  compartments  in  Figure  2.1  are  involved  with 
the  metabolic  processes.  The  metabolism  begins  when  the  hexane 
reaches  the  liver.  Much  of  the  hexane  is  oxidized  into  methyl- 
n-butylke tone .  The  metbyl-n-butylketone  is  then  distributed  to 
the  body  water,  metabolized  into  2,5-hexanedione,  or  metabolized 
into  other  substances  that  are  currently  not  of  interest.  Subsequently, 
the  2,5-hexanedione  in  the  liver  is  then  distributed  to  the  body 
water.  Finally,  the  2,5-hexanedione  in  the  body  water  is  then 
excreted  to  the  kidneys. 

Formulation  of  the  Model  Be  nations.  With  the  compartmental  model 
in  Figure  2.1,  the  corresponding  mathematical  relationships  are 


derived.  The  equations  involved  are  often  referred  to  as  mass 
balance  equations.  The  formulation  of  the  equations  closely 
follows  Mayberry’s  development,  although  Mayberry’s  formulation 
is  primarily  directed  towards  the  chemical  styrene  (MAXB82b:  CH  2) 
The  parameters  of  the  model  are  listed  and  defined  in  Table  2.1. 
Table  2.2  lists  either  the  numerical  values  or  the  defining 
relationships  used  for  the  parameters.  In  addition,  AFAMRL 
has  indicated  that  the  following  relationships  must  be  satisfied: 

QC=QF+QL  +  QS+QR 

CV  =  (QF*CVF  +  QL*CVL  +  QR'CVR  +  QS*CVS)/QC 

The  equations  of  the  compartmental  model  will  now  be  developed. 

Lung  Compartment.  For  the  lung  compartment,  the 
mass  balance  equation  is 

dAI  +  (amount  entering  from  the  venous  blood)  «  (2 

dAX  +  dAB  +  (amount  leaving  with  arterial  blood) 

This  equation  says  that  in  a  small  time  interval  dt, 

the  amount  of  hexane  gas  entering  the  lung  compartment 
from  inhaled  air  and  venous  blood  returning  from  the  body 
equals  the  amount  of  gas  absorbed  into  the  blood  stream 
plus  the  amount  of  gas  leaving  the  lung  compartment  in 
the  arterial  blood.  (MATB82b:  2.14) 


iP^agggasBssgies  |g§§j55e*6 


TABLE  2.1  Parameters  of  the  Open  Chamber  Hexane  Model 


s  amount  of  hexane  in  capillary  blood  in  the  lungs  (mg) 
s  amount  of  hexane  in  fat  (mg) 
s  amount  of  hexane  in  inhaled  air  (mg) 

=  amount  of  hexane  in  liver  (mg) 
s  amount  of  methyl-n-butylketone  in  liver  (mg) 
s  amount  of  2, 5-hexanedione  in  liver  (mg) 

3  amount  of  hexane  metabolized  into  me  thy 1-n-b  uty Ike  tone  (mg) 
s  amount  of  me  t hyl-n-buty Ike  tone  metabolized  into  2 , 5-hexanedione  (mg) 
s  amount  of  methyl-n-buty Ike tone  metabolized  into  other 
metabolites 

s  amount  of  hexane  in  rapidly  perfused  tissue  (mg) 
s  amount  of  hexane  in  slowly  perfused  tissue  (mg) 
s  amount  of  methyl-n-butylketone  in  body  water  (mg) 
s  amount  of  2 , 5-hexanedione  in  body  water  (mg) 
s  amount  of  hexane  in  exhaled  air  (mg) 

=  body  weight  of  all  mammals  in  the  te6t  chamber  (mg) 
s  concentration  of  hexane  in  arterial  blood  (mg/1) 

=  concentration  of  hexane  in  fat  (mg/1) 
a  concentration  of  hexane  in  inhaled  air  (mg/1) 

=  concentration  of  hexane  in  liver  (mg/1) 
s  concentration  of  methyl-n-butylketone  in  liver  (mg/1) 

3  concentration  of  2, 5-hexanedione  in  liver  (mg/1) 
s  concentration  of  hexane  in  rapidly  perfused  tissue  (mg/1) 
s  concentration  of  hexane  in  slowly  perfused  tissue  (mg/1) 
s  concentration  of  hexane  in  venous  blood  (mg/1) 

s  concentration  of  hexane  in  venous  blood  leaving  fat  tissue  (mg/l) 

s  concentration  of  hexane  in  venous  blood  leaving  liver  (mg/l) 

s  concentration  of  hexane  in  venous  blood  leaving  rapidly 

perfused  tissue  (mg/l) 

CVS  s  concentration  of  hexane  in  venous  blood  leaving  slowly 
perfused  tissue  (mg/l) 

CX  s  concentration  of  hexane  in  exhaled  air  (mg/l) 

K  s  rate  constant  for  the  excretion  of  2, 5-hexanedione  (1/hr) 

KFO  s  rate  constant  for  first  order  metabolism  (mg/l) 
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TABLE  2.1  (continued) 


KM1  =  enzyme  affinity  for  hexane  (mg/1) 

KM 2  =  enzyme  affinity  for  methyl-n-buty Ike tone  (mg/1) 

KM2I  =  constant  for  describing  interference  of  metyl-n-butylketone 
on  hexane  metabolism  (mg/1) 

KM3  =  enzyme  affinity  for  2, 5-hexanedione  (mg/1) 

LI  =  ratio  of  molecular  weights  of  methyl-n-butylketone  to 
hexane 

L2  =  ratio  of  molecular  weights  of  2 , 5-hexanedione  to  methyl- 
n-butylketone 

PB  =  blood  to  air  partition  coefficient 

PF  =  fat  to  blood  partition  coefficient 

PL1  =  hexane  liver  to  blood  partition  coefficient 

PL2  =  methyl-n-butylketone  liver  to  blood  partition  coefficient 

PL3  =  2, 5-hexanedione  liver  to  blood  partition  coefficient 

PR  -  rapidly  perfused  tissue  to  blood  partition  coefficient 

PS  =  hexane  slowly  perfused  tissue  to  blood  partition  coefficient 

QC  =  cardiac  output  (1/hr) 

QF  =  fat  compartment  blood  flow  rate  (l/hr) 

QL  =  liver  blood  flow  rate  (1/hr) 

QP  a  alveolar  ventilation  rate  (1/hr) 

QR  =  rapidly  perfused  tissue  flow  rate  (lAr) 

QS  a  slowly  perfused  tissue  flow  rate  (1/hr) 

RAMI  a  rate  that  hexane  metabolizes  into  methyl-n-butylketone  (mg/hr) 
RAM2  a  rate  that  methyl-n-butylketone  metabolizes  into  2,5-hexa-Mione 
(mg/hr) 

RAM3  =  rate  that  methyl-n-butylketone  metabolizes  into  other  substances 
(mg/hr) 

VF  a  volume  of  fat  compartment  (1) 

VL  a  volume  of  liver  compartment  (l) 

VMAX1  a  maximum  velocity  of  metabolism  for  hexane  (mgAr) 

VMAX2  a  maximum  velocity  of  metabolism  for  methyl-n-butylketone  (mg/hr) 
VMAX3  a  maximum  velocity  of  metabolism  for  2 , 5-hexanedione  (mg/hr) 

VR  a  volume  of  rapidly  perfused  tissue  compartment  (l) 

VS  a  volume  of  slowly  perfused  tissue  compartment  (l) 

W  a  volume  of  distribution  of  2, 5-hexanedione  (1) 


TABU2  2.2  Parameter  values  and  defining  relationships 


B V 

CA 

CF 

CL1 

CL2 

CL3 

CR 

cs 

cv 

CVF 

CVL 

CVR 

CVS 

CX 

KH1 

LI 

L2 

PB 

PF 

PL1 

PL2 

PL3 

PR 

PS 

QC 

QF 

QL 

QP 

OR 

QS 

RAMI 

RAM2 

RAM3 

VF 

VL 

VMAX1 

VR 

VS 

VW 


.3 

(QC*CV  +  CI*QP)/(QP/PB  +  QC) 

AF/VF 

AL1/VL 

AL2/VL 

AL3/VL 

AR/VR 

AS/VS 

(QF*CVF  +  QL*CVL  +  QR*CVR  +  QS*CVS)/QC 
AF/(VF*EF) 

AL1/(VL*PL1 ) 

AR/(VR*PR) 

AS/(VS*PS) 

CA/PB 

-36 

1.186 

1.14 

.8 

70.0 

2.3 

1.0 

1.0 

4.0 

1.26 

4.13  =  QF+QL  +  QS+QR 

.04*QC 

•3*QC 

3.85 

.48*QC 

•2*QC 

KFO*CVL*VL  +  VMAX1  *AL1/(VL*PL1  *KM1  *  ( 1  +  (AL2  +  AL3)/(VL*KM2I)) ) 
VMAX2*AL2/(AL2  +  VL*PL2*KM2  +  PL2*AL1/KM1 ) 
VMAX3*AL2/(VL*PL2*KM3  +  AL2) 

.04*  BW 
.04* BW 
1.8 
.05*BW 
.78*B W 
.12 


+  QL  +  QS  +  QR 
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With  the  following  relationships, 


dAI  *  CI*QP*dt 
dAX  S  CX*QP*dt 
dAB  a  dCA’dt 

(amount  entering  from  venous  blood)  =  QC*CV*dt 
(amount  leaving  with  arterial  blood)  =  AC*CA*dt 


equation  (2.1)  may  be  written  as 

CI*QP*dt  +  QC*CV*dt  =  CX*QP*dt  +  dAB  +  QC*CA*dt 

Rewriting  equation  (2.7),  we  have 

dAB  a  -QC*CA*dt  -  QP*CX*dt  +  CI*QP*dt  +  QC*CV*dt 

Differentiating  equation  (2.4)  with  respects  to  t  gives 

dAB/dt  a  dCA/dt 

Now  differentiating  equation  (2.8)  with  respects  to  t  and 
substituting  in  equation  (2.9)  gives 


dCA/dt  =  -QC*CA  -  QP*CX  +  CI*QP  +  QC*CV 

i 


(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 


(2.7) 


(2.8) 


(2.9) 


(2.10) 
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In  the  past,  researchers  at  AFAMRL  have  noted  that  the 
concentration  in  CA  reaches  equilibrium  rapidly  in  relation  to  the 
other  compartments.  This  observation  led  to  the  assumption  that 


dCA/dt  is  approximately  zero.  Hence,  equation  (2.8)  becomes 


QC*CA  +  QP*CX  =  CI*QP  +  QC*CV 


(2.10) 


Since  CX  =  CA/PB,  then  equation  (2.10)  is  written  as 


QC*CA  +  QP* CA/PB  =  Cl *QP  +  QC*CV 


(2.11) 


or,  after  rearranging. 


CA  =  (QC*CV  +  CI*QP)/(QP/PB  +  QC) 


(2.12) 


Fat.  Slowly  Perfused  Tissue,  and  Rapidly  Perfused  Tissue 
Compartments.  The  compartments  of  fat,  slowly  perfused  tissue 
and  rapidly  perfused  tissue  have  mass  balance  equations  of  a  similar 
form.  The  equation  for  fat  will  be  developed  first. 

For  fat,  we  have  the  following  mass  balance  equation. 


qp*CA*dt  =  QF*CVF*dt  +  dAF 


(2.13) 
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As  Mayberry  points  out,  equation  (2.1 3)  says  that  for  a  snail 
time  interval  dt,  the  amount  of  hexane  entering  the  compartment 
through  the  arterial  blood  is  equal  to  the  amount  of  hexane 
leaving  the  compartment  through  venous  blood  plus  the  amount  of 
hexane  absorbed  by  the  compartment. 

Equation  (2.13)  is  rewritten  as 

dAF/dt  =  QF*CA  -  QF*CVF  (2.1 

or,  since  CVF  =  AF/(VF*EF),  equation  (2.14)  becomes 

dAF/dt  a  QF*(CA  -  AF/(VF*FF))  (2.1 

Equations  (2.14)  and  (2.13)  a re  in  a  form  that  may  be  implemented 
on  the  computer  and  solved  using  a  numerical  integration  routine. 

The  differential  equations  are  similarly  derived  for  the  slowly 
perfused  tissue  and  rapidly  perfused  tissue.  For  the  slowly 
perfused  tissue,  the  equations  are 

dAS/dt  =  QS*(CA  -  CVS)  (2.1 

or 


dAS/dt  =  QS*(CA  -  AS/(VS*PS) ) 


(2.1 


is  '.-v*.  v*  *1 


For  the  rapidly  perfused  tissue,  the  equations  are 


dAS/dt  =  QR*(CA  -  CVR) 


(2.18) 


dAR/dt  =  QR* (CA  -  AR/(VR*PR)) 


(2.19) 


Hexane  in  Liver  Compartment.  The  mass  balance  equation 
for  hexane  gas  in  the  liver  is 


QL*CA*dt  =  dALI  +  dAMI  +  QL*AL1*dt/(VL*PLl) 


(2.20) 


where 


QIi*CA*dt  =  amount  of  hexane  entering  the  liver  compartment 
dALI  =  amount  of  hexane  absorbed  by  the  liver  tissue 

dAMI  s  amount  of  hexane  metabolized  into  methyl-n- 

butylketone 

QL*AL1*dt/(VL*PL1 )  =  amount  of  hexane  that  leaves  the 

liver  compartment  through  venous 
blood 


Equation  (2.20)  may  be  written  as 


dALI/dt  =  QL* (CA  -  CVL)  -  dAMI/dt 


(2.21) 


Scientists  at  AFAMRL,  however,  have  determined  that 
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dAM1/dt  =  RAMI  =  KPO*CVL*VL  + 


(2.22) 


VMAX1  *AL1/(VL*PL1  •KMI  *  (1  +  (AL2  +  AL3)/(VL*KM2I))) 


Substituting  RAMI  and  CVL  =  AL1/(VL*PL1)  into  equation  (2.21)  gives 


dALI/dt  =  QL*(CA  -  AL1/(VL*PL1 ) )  -  RAMI 


(2.23) 


Methyl-n-butylketone  in  Liver  Compartment.  The  mass 
balance  equation  for  methyl-n-butylketone  that  is  in  the  liver 


QL*AW1*dt/VW  +  dAM1*L1  =  dAL2  +  dAM2  +  dAM3  + 

QIi*AL2*dt/(VL*EL2) 


(2.24) 


where 


QL*AW1*dt/VW  =  amount  of  methyl-n-butylketone  entering 

the  liver  compartment  from  the  body  water 
dAM1*L1  =  amount  of  hexane  metabolized  into  methyl- 

n-butylketone  and  retained  in  liver 
dAL2  =  amount  of  methyl-n-butylketone  absorbed 

by  the  liver  tissue 

dAM2  s  amount  of  methyl-n-butylketone  metabolized 

into  2,5-hexanedione 

dAM3  =  amount  of  methyl-n-butylketone  metabolized 

into  other  substances 

QL*AL2*dt/( VL*PL2 )  =  amount  of  methyl-n-butylketone  that 

is  distributed  to  body  water 
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Writing  equation  (2.2*0  as  a  differential  equation,  we  hare 


-dAL2/dt  »  QL«AL2/(VL*PL2)  -  QL'JWI/W  -  (2.25) 

dAM1/dt*L1  ♦  dA.'Vdt  ♦  dAM3/dt 


But,  since 


dAMI/dt  a  RAMI 

dAM2/dt  a  RAM2  a  VMAX2*AL2/(AL2  ♦  VL*PL2*KM2  +  PL2*AL1/KM1) 
dAM3/dt  a  RAM3  =  VMAX3*AL2/(VL*PL2*KM3  ♦  AL2) 


we  may  write  equation  (2.25)  as 


-dAL2/dt  a  QL*AL2/(VL*PL2)  -  QL*AY1/W  -  RAM1*L1  +  (2.26) 

RAM2  +  RAM3 


2. 5-hexane di one  in  Liver  Compartment.  The  mass  balance 
equations  for  2,5-hexanedione  that  is  in  the  liver  compartment  is 


QL*AV2  *dt/VW  +  dAM2*L2  =  dAL3  +  QL*AL3*dt/(VL*PL3)  (2.27) 

where 

QL"AW2/VW  s  amount  of  2,5-hexanedione  entering  from 
the  body  water 

dAM2*L2  =  amount  of  2,5-hexanedione  metabolized  from 
methyl -n-buty Ike tone  and  retained  in  liver 
dAL3  a  amount  of  2,5-hexanedione  absorbed  by  liver  tissue 
QL*AL3/(VL*FI»3)  =  amount  of  2,5-hexanedione  distributed  to 
body  water 


Writing  equation  (2*27)  as  a  differential  equation,  substituting 
dAM2/dt  a  HAM2,  and  rewriting,  ve  hare 


dAL3/dt  a  QL*AV2/W  +  L2*RAM2  -  QL*AL3/(VL*PL3) 

Hethrl-n-butrlketone  in  Body  Water  Compartment.  The 
mass  balance  equation  for  the  amount  of  aethyl-n-butylketone 
in  body  water  is 


dAWI  =  QL*CL2*dt/(VL*PL2)  -  QJ,*AW1*dt/VW 

where 


AVI  a  amount  of  methyl-n-buty  Ike tone  in  body  water 
QL*CL2*dt/(VL*EL2)  a  amount  of  methyl -n-butylkstone 

entering  from  the  liver 

QL*AW1*dt/VW  a  amount  of  oethyl-n-butylketone  leaving 
the  body  water  and  entering  the  liver 


Writing  equation  (2*29)  as  a  differential  equation,  we  have 


dAV1/dt  =  QI,*CL2/(?r**L2)  -  QL-AW1/W 


?.s.haMnedione  in  Body  Water  Compartment.  The 
mass  balance  equation  for  2 , 5-hexanedione  in  body  water  is 


(2.28) 


(2.29) 


(2.30) 


dAW2  a  QL*C13*dt/(VL*PL3)  -  QL*AV2*dt/W  -  K*AV2*dt  (2.31) 


where 


JW2  =  amount  of  2 , 5-hexanedi one  in  body  water 
QL*CL3*dt/(VL*PL3)  =  amount  of  2, 5-hexanedi one  entering 

from  the  liver 

QL*AW2*dt/VW  =  amount  of  2 , 5-hexanedione  leaving  body 
water  and  entering  the  liver 
K*AW2*dt  s  amount  of  2,5-hexanedione  excreted  to  the 
kidneys 


Writing  equation  (2.31)  as  a  differential  equation,  we  have 


dAW2/dt  =  QL*CL3/(VL*PL3)  -  QL*AW2/VW  -  K*AW2 


(2.32 


Equation  (2.32)  marks  the  last  of  the  eight  differential 
equations  that  model  the  Open  Chamber  Hexane  Model.  Now,  we 
must  relate  these  differential  equations  to  the  measurements 
taken  during  the  experiment.  The  next  section  describes  the 
models  of  the  measurements. 


Measurement  Models.  During  the  experiment,  researchers 
have  collected  measurements  used  for  estimating  the  unknown 
parameters.  One  set  of  experimental  data  is  the  input  chamber 
concentration  level  of  hexane.  The  other  three  experimental 
data  sets  were  extracted  from  the  drug  concentration  levels 
inside  the  mammal.  The  observations  include  the  hexane  concentration 
level  in  the  blood,  the  concentration  level  of  methy-n-butylketone 


$ 


in  the  body  water,  and  the  concentration  level  of  2 , 5-hexane di one 
in  the  body  water.  The  three  measurement  models  give  the 
predicted  values  used  for  comparison  with  the  experimental 


observations. 

The  first  measurement  model,  concentration  of  hexane  in 
the  blood,  is  given  by 

7,  =  CA 

where  CA  is  defined  by  equation  (2.12).  The  second  measurement 
model,  the  concentration  level  of  methyl-n-butylketone  in  body 
water,  is  given  by 

y2  =  5*AW1 

Finally,  the  third  measurement  model,  the  concentration  level 
of  2,5-hexanedione  in  body  water,  is  given  by 

y3  =  5* AW 2 

iTwimown  Model  Parameters 


(2.33) 


(2.34) 


(2.35) 


The  values  for  seven  of  the  parameters  listed  in  Table  2.1  are 
unknown*  In  order  that  the  Hexane  Model  be  useful,  these  parameters  must 


be  estimated.  The  researchers  at  AFAMRL  have  requested  that, 
when  estimated,  the  parameters  satisfy  certain  constraints. 
These  constraints  ensure  that  the  estimates  will  make  physical 
sense.  Table  2.3  displays  the  unknown  parameters  and  their 
constraints. 


TABI£  2*3  Unknown  model  parameters 


2.4 

< 

KFO 

< 

4.4 

.2 

< 

VMAX2 

< 

2.0 

.1 

< 

KM2 

£ 

10.0 

1.0 

< 

VMAX3 

< 

5.0 

.1 

< 

KM3 

< 

10.0 

1.0 

s 

KM2I 

< 

1000.0 

.1 

< 

K 

< 

.6 

Experimental  Data 

The  data  was  obtained  from  AFAMRL.  Four  experiments 


were  run,  each  with  a  different  time  history  of  the  input 

concentration  level  of  hexane.  Each  experiment  was  started 

at  a  relative  time  point  of  zero.  The  data  is  given  in  Table  2.4, 


TABLE 

2.4(a) 

Experimental 

data: 

Run  #1 

Time  (hra) 

‘input 

* Hexane 

*MNBK  *2,5-HD 

6.0 

1.8 

1.3 

.7 

3.2 

6.5 

0.0 

— 

— 

3.6 

7.0 

0.0 

— 

— - 

3.6 

8.0 

0.0 

— _ 

— 

2.8 

10.0 

0.0 

... 

... 

2.0 

12.0 

0.0 

... 

... 

.8 

14.0 

0.0 

... 

... 

.4 

18.0 

0.0 

.2 

TABLE 

2.4(b) 

Experimental 

data: 

Run  #2 

Time  (hrs) 

Input 

Hexane 

MNBK  2,5-HD 

.5 

3.5 

1.6 

.8 

... 

1.0 

3.5 

1.4 

.9 

... 

2.0 

3.5 

1.7 

1.6 

— 

3.0 

3.5 

2.2 

1.8 

... 

4.0 

3.5 

1.6 

1.6 

... 

6.0 

3.5 

2.2 

2.4 

... 

6.5 

0.0 

— 

— 

5.6 

7.0 

0.0 

— 

— 

5-1 

8.0 

0.0 

— 

... 

6.1 

10.0 

0.0 

— 

... 

4.7 

12.0 

0.0 

— 

— 

4.2 

Note:  Input  =  input  chamber  concentration  level  (mg/l) 

Hexane  =  hexane  concentration  level  in  blood  (mg/L) 

MNBK  =  concentration  of  methyl-n-butylketone  in 
body  water  (mg/l) 

2,5-ND  =  concentration  of  2t5-hexanedione  in  body  water  (mg/l) 


TABLE 

2.4(c) 

Experimental 

data:  Run  if 3 

Time  (hrs) 

Input 

Hexane 

MNBK 

2,5-HD 

.5 

10.6 

1.6 

.8 

1.0 

10.6 

1.4 

.9 

— 

2.0 

10.6 

1.7 

1.6 

— 

3.0 

10.6 

2.2 

1.8 

— 

4.0 

10.6 

1.6 

1.6 

— 

6.0 

10.6 

2.2 

2.4 

6.5 

0.0 

— 

1.9 

5.6 

7.0 

0.0 

— 

— 

5.1 

8.0 

0.0 

— 

— 

6.1 

10.0 

0.0 

— 

... 

4.7 

12.0 

0.0 

— - 

- — 

4.2 

TABLE  2.4(d)  Experimental  data:  Run  #4 


5umnia^^_of__the_C108ed_Chamber_HexsLne  Model  Equations 
The  model  equations  are 

X1  a  dAF/dt  a  QF*(CA  -  CVF) 

X2  a  dALI/dt  »  QL*(CA  -  AL1/(VL*PL1 ) )  -  RAMI 
Xj  ■  dAL2/dt  =-QL*AL2/(VL*FL2)  +  QL*CW1/VW  + 

HAM1 *L1  -  RAM2  -  RAM3 

• 

X4  =  dAL3/dt  =  QL*AW2/VW  +  L2*RAM2  -  QL*AL3/(VL*PL3) 
*5  =  dAW1/dt  =  QL*CL2/(VL*PL2)  -  QL*AW1/VW 
Xg  =  dAW2/dt  =  QL*CL3/(V1*PL3)  -  QL*AW2/VW  -  K*AW2 
i,  =  dAS/dt  a  QS*(CA  -  CVS) 

Xg  =  dAR/dt  =  qR*(CA  -  CVR) 


The  measurement  equations  are 

y1  a  CA  =  (QC*CV  +  CX*QP)/(QP/EB  +  QC) 
y2  =  5*AW1 
y?  =  5*AW2 
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III.  Background 


System  identification  has  long  been  an  area  of  interest  and 
application.  Techniques  of  system  identification  are  intended 
to  help  develop  models  that  describe  the  phenomena  of  interest. 
These  models  are  usually  mathematical  relationships  among 
important  variables.  Typically,  these  models  describe  physical 
processes,  although  equations  of  more  abstract  processes  are  also 
popular.  The  models  may  range  anywhere  from  a  fully  defined  set 
of  equations  with  perhaps  a  few  unknown  parameters  to  a  model 
that  is  completely  unspecified.  However,  the  specification  of 
the  model  and  its  unknowns  will  determine  the  results  from 
estimation.  Thus,  consideration  must  be  given  to  modelling  the 
desired  phenomena  appropriately. 

Another  critical  issue  in  system  identification  is  that  of 
deciding  how  to  evaluate  the  agreement  between  the  theoretical 
model  and  the  observations.  This  measure  of  agreement  will  guide 
the  selection  of  values  for  the  unknown  parameters.  The  values 
that  give  the  optimal  measure  of  agreement  will  be  chosen  as 
the  estimates.  Methods  of  selecting  the  estimates  are  usually 


called  optimization  methods.  Researchers  must  also  be  concerned 
with  the  collection  of  experimental  data.  Sometimes,  control 


may  be  exercised  over  the  phenomena  of  interest.  In  that  case, 
researchers  nay  be  able  to  design  experiments  and  collect  data 
in  a  manner  that  will  provide  very  precise  estimates. 

This  chapter  presents  background  material  that  may  be  useful 
in  the  identification  of  parameters  of  the  Open  Chamber  Hexane 
Model.  The  topics  presented  are  modeling,  optimization  and 
experimentation.  Occasionally,  some  of  the  topics  will  be  mentioned 
briefly.  The  reason  is  that  each  of  the  topics  and  many  of 
their  subtopics  involve  concepts  that  are  very  detailed  and 
require  much  study  to  understand  fully.  Still,  these 
topics  are  mentioned  since  they  may  prove  useful  in 
indicating  where  we  might  go  to  pursue  further  an  important 
subject  for  future  applications. 

Stochastic  Modeling 

The  Open  Chamber  Hexane  Model  presented  in  Chapter  II  is 
a  deterministic  model.  The  uncertainties  that  may  be  involved 
in  the  model  have  not  been  addressed.  The  deterministic  model 
may  be  appropriate  if  the  uncertainties  are  negligible.  However, 
researchers  at  AFAMRL  have  indicated  that  modeling  these  uncertainties 
may  be  appropriate  if  at  all  possible  in  future  applications. 


Hence,  a  review  of  some  techniques  dealing  with  these  uncertainties 
may  be  beneficial. 

In  the  Hexane  Model,  uncertainties  derive  from  three  sources. 

The  first  source  stems  from  the  fact  that  not  sill  is  known  about 
the  true  metabolic  and  distribution  processes  of  drugs  in  the  mammal. 

In  addition,  not  all  of  the  processes  involved  are  being  modelled. 

Noise  variables  may  be  used  to  account  for  the  underlying  processes 
not  included  in  the  model.  The  second  source  of  uncertainty 
stems  from  the  fact  that  measurement  data  itself  is  corrupted 
by  the  measuring  equipment  and  procedures  followed  during  the 
experiment.  Finally,  the  results  are  expected  to  vary  since 
each  mammal  of  the  same  species  does  not  possess  the  same  physiological 
characteristics.  Estimating  parameters  of  dynamic  models  in  the 
presence  of  noise  requires  knowledge  about  statistics  and  stochastic 
processes.  This  section  reviews  some  of  these  concepts  for  both 
linear  and  nonlinear  models. 

Linear  Dynamic  Stochastic  Models.  A  linear  dynamic  system 
with  time-invariant  parameters  and  no  noise  may  be  written  as 

X(t)  =  F*X(t)  +b*U(t)  (3.1a) 

Z(t)  =  h*X(t)  (3.1b) 

For  the  case  of  linear,  time-invariant  pharmacokinetic  models, 


X(t)  represents  the  drug  concentration  levels  in  each  compartment 
U(t)  is  the  hexane  concentration  in  the  air  and  Z(t)  represents 
the  measurements  of  the  model.  F,  b,  and  h  are  constant  matrices 
that  describe  the  relationship  between  X(t),  U(t),  and  Z(t). 

One  way  to  account  for  model  uncertainties  is  to  add  noise 
terms  to  the  noise  free  model.  Adding  noise  terms  to  the 
equations  of  (3.1)  gives 

X(t)  =  F*X(t)  +  b*U(t)  +G*n1(t)  ( 

Z(t)  =  h*X( t )  +  n2(t)  ( 

In  general,  for  any  noise  terms  n^Ct)  and  n2(t),  the  model 
may  not  be  directly  solvable.  However,  if  n^(t)  and  n2(t) 
are  assumed  to  be  white  gaussian  noise  processes,  then  solutions 
for  X(t)  and  Z(t)  are  obtained  in  a  relatively  straight  forward 
manner. 

To  fully  characterize  n^t)  as  white  noise,  we  first  note 
that  it  is  a  continuous  function  of  time.  For  this  case,  n^ (t) 
has  infinite  variance  and  yet  it  is  uncorrelated  in  time.  It 
has  a  mean  of  zero  and  a  constant  power  spectrum.  Since  its 
power  spectrum  is  constant,  white  noise  is  typically  specified 
by  this  constant,  say  £.  Thus,  n^Ct)  may  be  specified  completely 
as  a  white  noise  process  with  a  power  spectrum  of  strength  £. 
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where 


i,(t)*dt  =  Br(t) 


Br(t)  a  brownian  motion  process 


The  noise  term  n2(t)  differs  from  n^t)  in  that  it  is  a 
discrete  function  of  time.  In  this  case  n2(t)  is  completely 
specified  by  a  mean  of  zero  and  a  covariance  matrix  R . 

The  solution  to  the  system  of  stochastic  differential 
equations  (3.2a)  becomes 

X(t)  =  W(t,t0)*X(t0)  +  J*  W(t,s)*b*U(s)*d8  + 


(3.3) 


(3.4) 


J  W(t,s)*G*dBr(s) 


where 

W(t,tQ)  =  solution  to  the  differential  equation 
d(W(t,tQ))/dt  =  F*W(t,tQ) 
and  £<t0,t0)  =  I 

The  solution  X(t)  is  a  gaussian  process  completely  characterized 
by  its  mean  and  covariance  matrix 
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mean  =  X(t)  =  W(t,tQ)*X(t0)  +  J  W(t,s)*b*U(s)*ds  (3.5a 

t0 

covariance  T 

matrix  =  P(t)  =  V(t,tQ)*P(t0)*W(t,t0)  ♦  (3«5b 

-t  T  T 

r  V(t,s)*G*£*G  *W(t,s)  *ds 

To  obtain  the  8olution  (3.4),  all  of  the  values  in  the  matricea 
F,  b,  G,  and  £  must  be  known.  However,  in  many  applications 
not  all  of  the  values  are  known.  In  aueh  cases,  parameter 
estimation  techniques  are  used  to  estimate  these  values. 

In  the  simple  case  with  unknown  values  only  in  F,  Eykhof f 
mentions  that  parameter  estimation  is  still  a  nonlinear  problem. 

"This  implies  that  all  approaches  to  this  problem  will  have  the 
characteristics  of  the  model-ad justment  technique,  i.e.  they 
will  be  of  an  iterative  type.”  (ETKH74:  446-7)  So,  even  for 
this  simple  case,  no  closed  form  solution  can  be  found  generally 
for  estimating  the  unknown  parameters. 

Even  though  closed  form  solutions  for  the  parameter 
estimates  may  not  be  available,  solutions  may  be  found  with  numerical 
techniques.  Unfortunately,  not  much  has  been  said  in  the  literature 


about  the  exact  probability  distributions  of  the  estimates. 

This  makes  it  difficult  to  determine  how  precise  our  estimates 
are.  However,  Goodwin  and  Payne  provide  an  equation  that  may 
be  used  to  estimate  the  precision  of  the  estimates.  (GOOD77:  99-102) 
The  Cramer-Rao  lower  bound  is  another  approximation  to  the 
precision  of  the  estimates.  (EYKH74:  27) 

Nonlinear  Dynamic  Stochastic  Models.  A  typical  representation 
of  nonlinear  stochastic  differential  equations  is  given  by 

X  =  F(X(t),U(t),£)  ♦  0*n(t)  (3.6) 

In  general,  even  when  the  unknown  parameters  q  are  assumed  known, 
the  solution  X(t)  to  this  equation  is  still  extremely  difficult 
to  find.  Usually  many  approximations  must  be  made  to  the  equation 
in  order  to  solve  it.  Springer  reports,  however,  that  the  theory 
of  H-functions  may  be  useful  in  dealing  with  such  nonlinear 
estimation  problems.  (SFRI79:  8)  Unfortunately,  it  was  not  possible 
to  investigated  further  the  use  of  H-functions  for  this  thesis  effort. 

A  technique  that  may  be  useful  in  estimating  the  parameters  of 
nonlinear  as  well  as  linear  models  is  the  jacknife.  A  version  of  the 
jacknife  was  developed  by  Quenouille  (QUEN56)  and  then  later  developed 
by  Takey  (TUKE58).  Mostellor  and  Tukey  state  that 


The  name  "jacknife"  is  intended  to  suggest  the  broad 
usefulness  of  a  technique  as  a  substitute  for  specialized 
tools  that  nay  not  be  available ,  just  as  the  Boy  Scout's 
trusty  tool  serves  so  variedly.  The  jacknife  offers 
ways  to  set  sensible  confidence  limits  in  complex 
situations.  (HQST77:  133) 


Because  of  its  wide  generality  in  estimating  parameters  and  their 
variances,  the  jacknife  technique  may  prove  to  be  a  valuable  tool 
in  system  identification.  However,  because  of  time  constraints, 
it  was  not  possible  to  apply  the  jacknife  to  the  Hexane  Model. 

Optimization 

Numerous  techniques  have  been  developed  for  numerical 
optimization.  Several  of  these  techniques  are  discussed  by 
Bard  (BAHD74)  and  fiykhoff  (£1X37*0.  Other  techniques  may  be 
found  by  searching  through  the  many  texts  and  articles  dealing 
with  optimisation  and  mathematical  programming.  Only  a  few 
of  these  techniques  will  be  mentioned  in  this  section. 

Most  optimization  methods  useful  in  parameter  estimation 
fall  under  the  categories  of  local  and  global  optimization. 

Also,  optimization  techniques  exist  for  dealing  with  the  added 
complexity  of  constraints. 


Parameter  estimation  of  pharmacokinetic  models  requires 


the  techniques  to  handle  global  optimization  with  constraints.  However 
global  optimization  techniques  and  methods  to  handle  constraints 
both  typically  build  upon  the  techniques  of  local  optimization. 
Therefore,  it  is  fitting  to  examine  first  some  of  the  literature 
dealing  with  unconstrained  local  optimization. 

Unconstrained  Local  Optimization.  Of  all  the  numerical 
optimization  methods  in  existence,  those  dealing  with  unconstrained 
local  optimization  are  probably  the  most  numerous  and  widely 
discussed.  However,  only  a  few  of  the  unconstrained  local 
optimization  techniques  will  be  discussed.  In  order  to  discuss 
these  techniques,  the  following  notations  will  be  used: 

£  »  the  vector  of  unknown  parameters 

E  a  the  error  function 

The  objective  of  parameter  estimation  is  to  minimize  the  error 
function  E  where  E  is  a  measure  of  the  distance  between  the 
observations  and  the  model's  predicted  values. 

Steenest  Descent  Method.  The  Steepest  Descent  algorithm 
is  given  by  (ETKH74:  158) 


£<i+l)  a  £<i)  -  J*jB/j<J 


t 


J>0 


(3.7) 


This  algorithm  updates  the  estimate  of  £(i)  each  iteration  with 
a  correction  factor,  ^  For  sufficiently  small 

J,  this  updated  value  will  reduce  the  value  of  S.  The  gradient 
SE/>£  indicates  the  direction  to  change  cj.(i)  in  order  to  reduce 
E  by  a  maximum  amount.  Unfortunately,  this  maximum  reduction 
of  E  may  occur  only  for  very  small  excursions  away  from  <j(i). 
Without  further  knowledge  about  the  function  E,  we  cannot  be 
sure  that  a  large  change  of  £(i)  in  the  direction  aE/5£ 
decreases  E  at  all. 

The  constant  J  should  be  chosen  so  that  the  correction  will 
yield  a  decrease  in  E  for  each  iteration.  If  J  is  too  large,  it 
is  possible  that  the  algorithm  trill  overshoot  the  minimum 
value.  On  the  other  hand,  if  J  is  too  small,  the  algorithm 
may  move  too  slowly  towards  the  minimum  and  waste  computer 
resources.  To  choose  the  best  value  for  J,  we  may  need  to  try 
many  values  and  choose  the  value  that  does  well  in  reducing 
the  error. 

The  method  of  Steepest  Descent  does  have  the  advantage 
of  guaranteed  convergence  to  a  local  minimum  (ETKH74:  161). 

This  is  not  necessarily  true  for  other  methods.  The  disadvantage 
of  the  Steepest  Descent  method  is  that  it  say  not  converge  as 
rapidly  as  these  other  methods. 


This  algorithm  is  given  by 


a<i+i)  =  a<i)  -  £=a(i) W-a£  £=£(i) 


(3.8) 


This  algorithm  can  be  related  to  the  Steepest  Descent  method 

r ,  ,3-1 


by  letting  J  *  *1™  HUB  auuea  imgniiBuoa 

of  second  deri retires,  the  Newton-Raphson  method  generally  updates 
q(i)  at  each  iteration  with  greater  reduction  in  the  error  function 
E.  But,  to  ensure  such  greater  reduction,  the  Newton-Raphson 
method  should  hare  a  starting  value  close  to  the  minimum. 

As  noted  by  Eykhoff ,  "evaluation  of  the  second  derivative 
is  quite  a  troublsome  job."  (EIKH?4:  161)  However,  if  the  error 
function  E  is  a  weighted  least  squares  error  function,  then 
the  simpler  Gauss-Newton  method  may  be  used  to  approximate 


..  Vith  this  added  information 


the  Newton-Raphson  method. 


L.  For  a  weighted  least  squares 


error  function,  the  Newton-Raphson  method  cam  be  approximated 
by  the  Gauss-Newton  method,  fykhoff  provides  the  derivation 
for  the  case  of  continuous  measurements.  For  this  case,  the  error 
function  is 


r. 


<t)  •2*®<t)*dt 


(3.9) 


where 


e.  =  0(t)  -  i(£,t) 

£  a  a  positive-definite,  symmetric  weighting 

0(t)  *  measurement  data 

l<a,t)  a  the  model's  predicted  values 


From  the  error  function,  we  have  the  following: 


T 

3E/a£  a  -  f  ajr(£,t)  /aq*£*e*dt 


(3.10) 


a^E/sgp^  a  f ai(a,,t)T/aa*£*^l(a.,t)/a£T*dt 

V  A 


(3.11) 


Substituting  equations  (3.10)  and  (3.11)  into  (3.8)  gives 
the  following  Gauss-Newton  algorithm  (EIXH74:  161) 


£(i+l)  a  £(i 


o  4f 

LJc 


(32,T /«£)  *£*  (aj/9^  )*dt]  * 


&=a(i) 


(3.12) 


Oi  /aa)*£*e*<it  m(1) 


The  advantage  of  this  algorithm  over  the  Nevton-fiaphson 
algorithm  is  that  no  second  derivatives  need  to  be  calculated. 
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Hence,  the  computations  are  simplified.  Unlike  the  method  of 
Steepest  Descent,  however,  this  method  does  not  guarantee  convergence 
to  a  local  minimum.  On  the  other  hand,  when  close  to  the  minimum, 
the  Gauss-Newton  method  converges  faster  than  the  Steepest 
Descent  method. 

Box-Kanemaau  Interpolation  Method.  A  particular  algorithm 
described  by  Beck  and  Arnold  is  the  Box-Kanemasu  algorithm. 

The  Box-Kanemaau  algorithm  is  simply  a  modification  of  the 
Gauss-Newton  method.  ( BECK 77 i  362)  Beck  and  Arnold  report 
that  the  Box-Kanemaau  algorithm  is  superior  to  most  other 
optimization  techniques.  They  also  point  out  that 

Bard  [Bard,  p.  Ill]  appears  to  favor  a  modification 
of  the  Box-Kanemaau  method  which  he  calls  the 
interpolation-extrapolation  method.  (HECK77:  373) 

• 

Furthermore,  they  state  that  the 

Box-Kanemaau  method  is  recommended  in  order  to  be  more 
certain  of  finding  the  parameter  values  minimizing  S 
the  error  function,  even  though  in  some  cases  the  Gauss 
method  would  be  more  efficient.  (BECK 77:  378) 

Beck  and  Arnold  conclude  that  the  Gauss-Newton  method  and  the 
Box-Kanemasu  method  both  work  well  in  parameter  estimation 


problems. 


Conjugate  Gradient  Method.  The  Conjugate  Gradient 


method  is  presented  by  Fletcher  and  Beeves  (FIET64).  Chattergy 
and  Vismer  report  that  the  Conjugate  Gradient  method  converges 
faster  than  the  Steepest  Descent  method  near  the  optimum  point 
(CHAX78:  332)*  However,  Maybeck  points  out  that  "the  lit  arature 
reports  certain  undesirable  numerical  characteristics  of  this 
algorithm  for  the  parameter  estimation  application."  (MAYB82a:  82) 

It  appears  that  perhaps  not  much  can  be  said  about  the  performance 
of  the  Conjugate  Gradient  method  unless  we  implement  it  specifically 
on  the  model  of  interest  aand  evaluate  its  performance. 

Derivative-Free  Methods.  The  optimization  techniques 
mentioned  so  far  required  that  at  least  the  derivative  of  the 
error  function  be  calculated.  Sometimes,  however,  it  may  be 
too  difficult  to  find  analytically  tbs  derivatives.  In  this 
case,  a  derivative-free  method  may  prove  useful.  Two  derivative- 
free  techniques  described  in  this  section  are  the  method  of 
finite  differences  and  the  method  of  direct  search. 

With  the  method  of  finite  differences,  the  algorithms  discussed 
previously  may  still  be  used.  The  only  modification  involves 
approximating  the  derivatives  in  these  algorithms.  For  example, 
all  of  the  algorithms  required  that  the  gradient  be  known. 

A  finite  difference  approximation  to  this  vector  is  composed  of 
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elements  found  by  (BARD74:  117) 


ai/aqi  =  a£/a  (3.13) 

where 

=[E(q1*q2*  *•*  *  qi*h*  »qn*  “  (3.14) 

4 

and  is  sufficiently  small 

Similar  approximations  can  also  be  made  for  the  second  derivative. 

Other  derivative-free  techniques  are  the  direct  search 
techniques.  These  techniques  include  the  first  and  second 
Posell  methods,  Hooke  and  Jeeves  pattern  search,  Rosenbrock's 
method,  and  the  simplex  method  of  Holder  and  Mead,  to  name 
a  few.  (BAHD74:  120) 

Direct  search  methods  have  performed  well  in  some  cases, 
but  Bard  notes  their  inadequacy  for  parameter  estimation. 

Specifically,  Bard  states 

Our  own  experience,  however,  has  been  disappointing; 
gradient  methods,  even  using  finite  difference  approximations 
have  outperformed  direct  search  methods  on  all  but  the 
most  trivial  parameter  estimation  problems,  both  in  reliability 
and  speed  of  convergence.  (BARD74:  119-120) 


In  addition,  Banks  states  that 


These  methods,  which  make  use  of  only  functional 
evaluations,  are  attractive  in  some  cases  if  one 
suspects  that  the  function  J  to  be  minimized  is 
not  smooth,  but  they  are  slow  and  usually  quite 
inefficient  when  highly  accurate  solutions  are 
desired.  Furthermore,  they  have  been  developed 
heuristically  and  no  proofs  of  convergence  have 
been  given.  (BANK80:  13-14) 

In  light  of  these  comments,  no  attempt  is  made  to  explain 
in  detail  the  direct  search  algorithms.  The  reader  may  consult 
the  bibliography  given  by  Bard  for  further  references  to  these 
methods. 

Unconstrained  Global  Optimization.  The  optimization  methods 
discussed  up  to  this  point  seek  only  to  find  the  local  optimum 
of  the  objective  function.  As  a  result,  the  estimate  obtained 
by  these  methods  may  be  a  local  minima  and  not  global.  In  an 
attempt  to  correct  these  deficiencies,  global  optimization 
methods  have  been  developed.  No  known  technique  can  guarantee 
a  globally  optimal  solution.  However,  these  techniques  will 
increase  the  chances  of  obtaining  a  global  optimum. 

Bard  explains  that  most  of  the  local  minimization  techniques 
will  reach  a  global  minimum  provided  they  have  appropriate  starting 
values.  (BAHD74:  120)  Therefore,  a  scheme  for  selecting  good 
starting  values  should  increase  the  chance  of  finding  the  global 


minimum.  Two  such  schemes  explained  by  Bard  are  the  grid  search 
and  the  random  search  methods. 

In  the  grid  search  method,  a  feasible  region  is  first 
defined.  This  usually  requires  aprior  knowledge  about  the 
model  parameters.  Typically,  this  feasible  region  is  represented 
by  upper  and  lower  bounds  of  the  unknown  parameters.  The 
feasible  region  then  becomes  rectangular.  We  conduct  a  grid 
search  by  selecting  many  points  in  the  feasible  region  and 
evaluating  the  error  function  at  each  of  the  points.  The  point 
that  has  the  smallest  error  function  value  is  then  used  as 
the  starting  value.  (BABD74:  121) 

In  the  random  search,  test  points  are  selected  randomly. 

Then  the  point  that  has  the  smallest  value  in  the  objective 
function  is  used  as  the  starting  value.  (BAKD74:  121) 

Variations  of  both  the  grid  search  and  random  search  may 
be  developed.  For  example,  instead  of  using  only  one  point  for 
the  starting  value,  we  may  use  three  points.  Then  an  optimization 
algorithm  may  be  applied  with  each  of  these  values  as  starting 
values.  Clearly,  if  all  three  solutions  are  equal,  then  we 
may  be  confident  that  the  estimate  is  a  global  minimum.  However, 
if  the  three  solutions  differ,  then  this  indicates  that  many 
local  mini nrums  exist. 

Another  algorithm  given  by  Bremermann  also  attempts  to 


obtain  the  global  optimum.  (BREM70)  This  method  uses  random 
search  schemes  in  conjunction  with  interpolation  of  the  objective 
function  with  a  fourth  degree  polynomial.  Lehman  Stark  report 
that  with  twelve  test  functions,  the  Bremermann  algorithm  was 


Faster  than  fourteen  other  algorithms  when  used  on 
parameter  spaces  of  dimension  greater  than  four  ...  . 
Additionally,  the  random  search  has  the  advantage 
of  global  optimization:  the  program  does  not  get 
stuck  at  local  extrema,  as  gradient-search  techniques 
are  prone  to  do.  (LEHM81:  117-118) 


We  should  note  that  this  statement  errs  in  implying  that  this 
algorithm  finds  the  global  optimum.  Rather,  the  algorithm 
at  most  only  increases  the  chance  of  finding  the  global  optimum. 

Constrained  Optimization.  As  indicated  by  the  researchers 
of  AFAMRL,  the  unknown  parameters  of  the  Hexane  Model  must 
satisfy  certain  inequality  constraints.  Thus,  methods  must  be 
explored  that  ensure  our  estimates  will  satisfy  those  constraints 
Chattergy  and  Wismer  explain  that  (CHAT78:  163) 


There  are  many  algorithms  which  have  been  proposed  to  solve 
constrained  optimization  problems  using  gradient  techniques. 
Although  these  algorithms  are  too  numerous  to  attempt  a 
comprehensive  treatment  here,  they  can  generally  be  divided 
into  two  categories: 


(1)  boundary-following  methods  and 

(2)  penalty-function  methods. 


Under  the  category  of  boundary  following  methods  are  the 
method  of  feasible  directions  and  the  gradient -projection  method. 

For  a  further  discussion  of  these  methods,  the  reader  may  refer 
to  Rosen  (RQSE60;  R0SE61)  and  Zoutendijk  (Z0UT60). 

The  penalty-function  methods  for  handling  constraints  are 
relatively  simple.  Penalty  functions  are  used  "to  convert 
constrained  optimization  problems  into  unconstrained  problems. 

Once  the  unconstrained  problem  is  formulated,  direct  methods 
for  solving  this  type  of  problem  can  be  employed."  (CHAT78:  197) 
Thus,  the  penalty-function  methods  have  the  advantage  that  the 
unconstrained  optimization  methods  may  be  applied  directly 
without  modification. 

Experimental  Design  for  Dynamical  Systems 

Goodwin  and  Payne  state  that  "For  dynamic  systems,  experiment 
design  includes  choice  of  input  and  measurement  points,  test 
signals,  sampling  instants  and  presampling  filters."  (MEHR76:  252) 
These  topics  each  appear  to  be  relatively  complex.  A  detailed 
investigation  of  experimental  design  procedures  was  not  possible 
within  the  time  constraints  of  this  research.  Instead,  this  section 
will  provide  a  guide  to  some  of  the  literature  on  experimental 


design  for  dynamical  systems.  The  two  topics  that  will  be  discussed 
in  this  section  are  (1)  input  design  and  (2)  choice  of  sampling 


instants 


Siaha  and  Kuszta  describe  three  methods  for  the  optimal 
design  of  inputs.  They  are  (SINH83 :  230) 

a)  method  based  on  properties  of  Fisher's 
information  matrix 

b)  method  based  on  system  sensitivity 
function  optimization 

c)  method  based  on  Gagliardi  theorem  - 
the  "direct"  method 

In  addition  to  these  methods  are  input  designs  that  aid  in 
discriminating  among  many  possible  models.  Together,  these 
topics  form  a  large  part  of  the  methods  for  input  design. 

The  method  involving  Fisher's  information  matrix  theory  is 
sometimes  referred  to  as  the  frequency  domain  design  of  input 
signals.  (SINH83:  231)  This  method  is  also  explained  by  Mehra 
(MEHR76:  211-229)  and  Goodwin  and  Payne  (GOOD77:  133-157). 

Sinha  and  Kuszta  note  that  Fisher's  information  matrix  approach 
pertains  to  linear  systems  only. 

The  method  using  sensitivity  functions  may  also  be  used 
to  design  the  inputs  to  the  model.  Sensitivity  functions  are 
defined  aa  X.  ^  3  a^/aq^  where  n^  is  a  state  variable  and  q.. 
is  an  unknown  parameter.  Beck  and  Arnold  explain  how  the 


sensitivity  functions  affect  the  convergence  of  the  Gauss-Newton 
method.  In  particular,  the  matrix  of  sensitivity  functions  should 
be  invertible  near  the  region  of  the  minimum.  Should  the 
determinant  of  the  sensitivity  matrix  become  small,  then  "The 
minimum  point  will  not  be  very  pronounced.”  (BECK 77*  347)  In 
this  case,  we  might  expect 

slow  convergence  of  the  Gauss  method.  For  this  reason 
it  is  important  to  examine  the  sensitivity  coefficients 
over  the  region  of  interest.  (BECK77*  348) 

Furthermore, 

For  effective  nonlinear  estimation,  the  careful 
examination  of  these  sensitivity  coefficients  is 
imperative.  (HECK 7? :  349) 

A  method  of  using  sensitivity  functions  to  design  inputs 
is  explained  by  Lehman  and  Stark.  They  recommend  that  normalized 
plots  of  the  sensitivity  functions  be  compared.  From  these  plots, 
we  know  which  state  variables  are  affected  the  most  by  the 
unknown  parameters  at  any  time.  (LEHM81:  113-117)  Hence,  the 
most  precise  estimates  may  be  gained  by  gathering  data  when  the 
sensitivities  are  the  largest.  Beck  and  Arnold  describe  how  we 
might  normalize  the  sensitivity  functions.  First,  let  be 


the  sensitivity  function  of  the  parameter  q .  at  time  i.  Also, 

3 

n^  represents  state  variables  and  tf^(i)  is  the  standard  deviation 
of  the  measurement  at  time  i. 


For  single  response  cases  vith  approximately  constant 
standard  deviations  of  the  measurements  it  is  convenient 
to  examine 

X.  .  =  q  .*X.  .  =  q.*9n./3q. 

13  3  13  3  1  h3 
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Note  that  X.  .  has  the  units  of  n.  Then  the  magnitude 
of  each  sensitivity  can  be  compared  with  the  others  as 
well  as  with  n  itself. 

For  multiresponse  cases  it  is  often  more  meaningful 
to  plot 


Xj^(i)  =  (q^/&-k(i))*(3nk(i)/^q;j) 

which  is  dimensionless.  (Note  i  refers  to  '’time1',  j  to 
parameter,  and  k  to  response.) 


Choosing  inputs  to  maximize  the  sensitivities  may  be  considered 
as  a  mathematical  optimization  problem.  Using  calculus  of 
variations  techniques,  the  problem  is  reduced  to  a  two-point 
boundary  value  problem.  Me hr a  (MEHR76:  211-249)  and  Sinha  and 
Kuszta  (SINH83:  235-237)  examine  the  methods  of  input  optimization 
for  linear  systems.  Kalaba  and  Spingarn  consider  the  case  for 
nonlinear  systems.  (KALA82:  225-379) 


Sinha  and  Kuszta  provide  a  third  approach  to  input  optimization 
known  as  nthe  direct  method  of  optimal  input  design. "  (SINH83:  237-258) 


This  method  is  developed  only  for  linear  time  invariant  systems. 
Goodvin  and  Payne  provide  a  good  development  on  input  optimization. 
(GOOD77:  124-157) 

The  input  optimization  topics  have  so  far 


Been  concerned  with  the  design  of  experiments  for 
accurate  parameter  estimation  within  a  model  of 
specified  structure.  However,  it  is  sometimes 
the  case  that  there  are  two  or  more  rival  models 
and  the  purpose  of  the  experiment  may  be  to 
determine  which,  if  any,  of  the  models  are 
adequate.  (GOOD77J  167) 


The  topic  of  experimental  design  for  model  discrimination  is 
considered  by  Bard  (BABD74:  266-269),  Beck  and  Arnold  (BECK 77:  464- 
473),  and  Goodwin  and  Payne  (GOOD77:  167-172). 

A  final  topic  under  experimental  design  considered  in  this 
chapter  is  the  choice  of  sampling  intervals.  Goodwin  and  Payne 
state  in  the  book  edited  by  Mehra  that 


for  those  problems  where  the  choice  of  sampling  instants 
and  pre-sampling  filters  does  arise,  it  has  been 
recognized  that  this  choice  often  has  a  significant  effect 
on  the  information  return  from  the  experiment  ...  . 

The  choice  is  particularly  critical  when  the  available 
computer  storage  and  analysis  time  are  limited. 

In  these  cases,  the  experiment  should  be  designed 
to  maximize  the  average  information  content  of  each 
sample.  (MEHR76:  252) 


The  choice  of  sampling  intervals  may  therefore  be  an  important 
consideration  for  pharmacokinetic  models  since  experiments 
are  time  limited.  Goodwin  and  Payne  discuss  some  methods  in 
choosing  sampling  intervals t  although  these  methods  seem  to 
apply  only  to  linear  systems.  (GOOD77:  157-16?;  MEELR76:  251-287) 


IV.  Computer  Algorithm  Development  for  Parameter  Estimation 


Some  of  the  parameter  estimation  methods  discussed  in 
Chapter  III  are  developed  in  this  chapter.  Specifically*  the  optimiza¬ 
tion  algorithms  which  were  implemented  on  the  computer  are  described. 
The  model  being  considered  is  the  deterministic  representation  of  the 
Open  Chamber  Hexane  Model  presented  in  Chapter  II.  Parameter 
estimation  of  deterministic  models  has  been  the  approach  taken 
by  AFAMHL  in  the  past  and  is  the  approach  they  currently  recommend 
now. 

To  develop  the  optimization  routines*  the  error  function  was 
first  defined.  Once  the  error  function  was  defined*  computer 
algorithms  were  then  designed  and  programmed  to  estimate  the 
parameters.  Next*  the  computer  algorithms  were  tested.  Once 
they  were  verified  to  work  correctly,  the  computer  programs 
were  then  compared. 

Weighted  Least  Squares  Criterion 

The  error  function  used  for  estimation  is  the  weighted 


1 


least  squares  function.  The  primary  consideration  in  this  choice 
is  the  approval  of  the  sponsors  themselves.  However,  in  addition 
to  meeting  AFAMRL's  approval,  this  error  function  has  another 
advantage.  It  does  not  require  any  statistical  assumptions. 

(EYKH74:  39)  This  is  important  since  the  noises  in  the  Hexane 
Model  are  not  considered  in  this  thesis  effort. 

Specifically,  the  weighted  least  squares  error  function 
used  for  parameter  estimation  is  given  by 

[(0.  (k)  -  y.  .(k)  )/0.  .(k)l  (4.1) 

i=1  j-i  k=1  1  J  J  J 

where 

i  =  index  over  the  four  experimental  mins 
j  =  index  over  the  three  measurements 
k  =  index  over  the  fifteen  observation  time  points 
0_(k)  *  experimental  observations 

y^tk)  a  predicted  values 

One  important  issue  in  evaluating  the  error  function  is 
that  of  finding  y^Ck).  As  shown  in  Chapter  II,  y^Ck) 
derived  from  the  solution  of  the  differential  equations  of  the 
Hexane  Model.  A  great  deal  of  Mayberry's  work  dealt  with  the 
solution  of  the  differential  equations  of  the  pharmacokinetic 
models.  A  major  problem  noted  by  Mayberry  is  that  the  pharmacokinetic 


according  to  the  equation 


a(i+D  =  a(i)  -  J*3E/»a|a=a(i) 


(4.2) 


With  a  step  size  J  specified,  the  algorithm  then  begins  from 
some  initial  value  a(0)  provided  by  the  user.  The  only  other 


value  required  for  each  iteration  is  sE/s 


One  approach  to  computing  VE/dg.  follows  by  first  defining 


ei;.(k,q)  =  (0i;.(k)  -  yi;j(k,a))/Oij.(k) 


Thus,  the  error  function  may  be  written  as 


(4.3) 


E(a)  =  (ei;j(k,a)) 


(4.4) 


i=  i  .i=i 


Differentiating  equation  (4.4)  with  respect  to  q^,  we  see  that 


2*Ci  *(aei.(k,a)/9qm)  (4.5) 

i=1  j=1  k=1  J  J 


Differentiating  equation  (4.3)  with  respects  to  q^,  we  see  that 
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3e.j(kta)/jqm  =  (-»yi;j(lc,a)/jqja)/Oi  (k) 


(4.6) 


Substituting  equation  (4.6)  into  equation  (4.5),  we  have 


(4.7  > 


Hence,  from  equation  4.7,  we  now  see  that  the  computation  of 
3E(c[)/?£  reduces  to  the  computation  of  3y\  and  e„. 

Before  evaluating  3yi;.(k,£)/aqm,  we  first  recall  that  y^k,^) 
is  defined  by  equations  (2.33),  (2.34),  and  (2.35)»  In  other  words, 
the  measurements  yi ^(k,^)  are  related  to  the  state  variables 
xnj(ki£)*  To  evaluate  ay^(k,£^«gt,  we  thus  need  to  evaluate 
«Jxnj(k,<j)/3£.  The  procedure  for  doing  this  is  to  first  start  with 
the  original  differential  equations  of  the  Hexane  Model.  This  is 
given  by 


xuj(k,a)  =  fn(X,k,a) 


(4.8) 


Nov,  differentiating  equation  (4.8)  with  respect  to  each  of 
the  unknown  parameters,  we  have 


a(Xn;.(k,£))/3qffl  =  JfgQC'k,^)/^  (^«9) 

But, 

aX^k^/aq^  =  d(axn^(k,£)/aqm)/dt  (4.10) 

Thus,  the  sensitivity  function  sX^tkj^Vaq^  may  be  expressed 
in  terms  of  the  differential  equation 

d(3Xn;.(k,£)/aqm)/dt  =  af^X^jO^/aq^  (4.11) 


I 


i 


The  significance  of  equation  (4.11)  is  that  the  solution 
to  it  is  the  sensitivity  function  which  can  be  evaluated  with 
the  same  numerical  integration  routine  used  to  evaluate  the 
differential  equations  of  the  Hexane  Model.  However,  a  problem 
with  equation  (4.11)  arises  when  we  consider  the  number  of 
equations  it  represents.  To  use  equation  (4.11',  we  must 
differentiate  each  of  the  8  equations  MX,*,^)  7  times 
the  7  unknown  parameters.  This  gives  56  differential  equations 
that  must  be  derived.  Had  fQ(X,k,£)  been  linear,  deriving 
the  56  differential  equatoons  may  not  have  been  difficult. 
However,  since  the  Hexane  Model  is  nonlinear,  the  derivation 
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of  the  56  differential  equations  is  very  tedious.  Thus,  another 
method  for  evaluating  "3E( was  used. 

The  values  for  9E  (_£)/?£  were  approximated  using  the  finite 
difference  method  presented  in  Chapter  III.  Although  this 
method  is  oqly  an  approximation ,  it  has  the  advantage  that  no 
additional  differential  equations  were  required.  Thus,  only 
the  original  model  equations  were  needed  to  evaluate 
9E(c[)/fc£. 

As  implemented  on  the  computer,  the  Steepest  Descent  method 
accepts  the  new  value  £(i+1 )  as  the  current  estimate  only  if 
it  produces  an  error  sum  of  squares  value  less  than  that  of 
<L(i).  If  ^(i+1)  is  rejected,  the  step  size  J  is  adjusted 
to  a  smaller  value  and  the  iteration  begins  once  again.  However, 
once  J  becomes  less  than  some  very  small  value,  the  iteration 
stops  since  we  are  close  to  the  local  minimum. 

To  satisfy  the  boundary  constraints  placed  on  the  parameters 
additional  logic  was  included  into  the  computer  program.  When 
outside  its  limits,  the  parameter  is  replaced  by  its  constraint 
value.  Although  this  method  seems  to  work  well  for  the  Hexane 
Model,  there  exist  other  modifications  to  the  Steepest  Descent 
algorithm  that  will  enable  it  to  deal  with  the  constraints. 

The  algorithm  implemented  on  the,  computer  is  described  in 
Figure  4.1. 


Steepest  Descent  Algorithm 


with  a  starting  point  £  f  evaluate  ssq,  the  sum  of  squares  error 

DO  WHILE  ((step  size  J  .GT.  .000001)  .AND.  (counter  .IE.  max_count) 

•AMD.  (at  least  one  improvement  in  the  last  3  iterations)) 

FOR  i=1  to  7  DO 

q_temporary  =  q^ 

q.  =  1.00001*q. 
nx  1 

evaluate  ssq_new  for  the  new  vector  £ 

=  q__temporary 

AE/Aq^  =  (ssq_new  -  ssq)*100000/q^ 
q__test^  =  qi  -  J*AE/Aqi 
IF  (q__test^  .LT.  lowerbound^)  THEN 
q_test^  =  lowerbound^ 

ELSE  IF  (q_testi  .GT.  upperboundj  THEN 
q__test^  =  upperbound^ 

ENDIF 
ENDFOR  i 

evaluate  ssq_test  for  the  new  estimates  q  test 
IF  (ssq^test  .LT.  ssq)  THEN 
FOR  i=1  to  7  DO 
q^  =  q_test^ 

ENDFOR  i 
ssq  =  ssq_test 

record  that  an  "improvement"  occurred  this  iteration 
ELSE 

record  that  "no  improvement"  occurred  this  iteration 
reduce  J  to  J/10 
ENDIF 

increment  counter 
ENDWHILE 

the  final  estimate  is  £  with  a  sum  of  squares  value  ssq 

FIGURE  4.1  Steepest  Descent  Algorithm 
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Gauss-Newton  Method.  To  derive  the  Gauss-Newton  algorithm  for 
the  Hexane  Model,  we  start  with  the  error  function 


K(cl)  = 


i=1  j=T  k= 


&2 

(e.(k,£)) 

A  J 


(4.4) 


Next,  a  first  order  Taylor  series  approximation  to  (e.  .(k,c[)) 

1  3 

is  substituted  into  equation  (4.4)  to  get 


E(a)  =  bbh\  e,-  .0c,£O)  + 

i=i  j=i  k=i  L  J 


(4.13) 


C  (3e  (k,£°)/3qm)*AqJ 

m=1  0  J 


where 


q°  =  some  nominal  value 


Now  we  minimize  the  error  function  by  takiuf  derivatives  with 
respect  to  £  and  setting  the  result  to  zero.  This  gives 


E<a>/  q  =  t2  fZ  £  2’f ♦  C4 

i=1  j=1  k=1  [  13 

*3ei;j(k,£0)/^qn 


(4.14) 


for  n=  1 |  2^  e««  |  7 


1  •  *  * -  ji  * .  *  •  *j*  *,■  m*  ,*  .*  .  *  «*  «  «*  ^  «"  .  * ,  •«  •  ’  •  *  %*  ^  \  "»  **  i  *•  *. 

*.  *-  .*1  *•  •  -  •  .  *.  V  V  -p  ^  a*  ^  -.V  -  «  •  “  ■  ‘  A  “ 


Dividing  equation  (4.14)  by  2  and  rewriting,  we  then  have 


T2  ~ei  1.(k«q°)*(3e.  .(k.q°)/aq  )  =  (4 

i=1  k=1  A  J  1  o  11 

C  ^(6  <8.lj(k,a°)/5qn)-  («1J(k,a°)/»qB)*AJm 


But,  equation  (4.15)  nay  be  rearranged  to  get 


Y2  ^2^2  ^2  ^e..(k,£°)/3q  )*(ae.,(k,£°)/aq  )*^q 


i=1  j=1  k=1  m=1 


Equation  (4.16)  is  now  in  a  fora  that  we  can  represent  using 
matrix  equations.  The  advantage' of .equation  (4.16)  is 
that  we  can  now  solve  for  <j^  explicitly. 

In  order  to  solve  for  q^.  we  first  define  the  following 
relationships 


Da=  LUC  ei  j(k»a  )*(* ®i j  (k,a° )/3<i  ) 
i=1  j=1  k=1  13  13 


0«ij(k,a°)/»qil)*(»«i;j(k,a0)^qiii)  (4.18) 

Thus,  in  matrix  form  equation  (4.16)  may  be  written  as 

=  £*  1  =  £*(=L  -  £°)  (4.19) 

Solving  for  £,  we  have 

1  =  1  -  (4.20) 

With  £  a  £(i+1)  and  £°  a  £(i) ,  the  Gauss-Newton  algorithm  for 
the  Hexane  Model  is  written  as 

a(i+i)  =  a(i)  -  c-1*D|^(i)  (4.2i) 

Equation  (4.21)  is  the  specific  algorithm  used  for  estimating 
the  parameters  of  the  Hexane  Model.  The  algorithm  starts  with 
an  initial  value  of  £(0).  The  partial  derivatives  could  have 
been  determined  directly  from  the  sensitivity  functions. 

But,  as  was  shown  for  the  Steepest  Descent  method,  it  is  very 
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tedious  to  derive  the  differential  equations  of  the  sensitivity 
functions.  Therefore ,  the  finite- difference  method  is  used 
to  approximate  the  partial  derivatives.  We  do  note  that  equation 
(4.21)  is  similar  to  equation (3. 12)  except  that  equation  (4.21) 
deals  with  discrete  measurements  instead  of  continuous  measurements. 

The  method  used  for  dealing  with  constraints  is  the  same 
as  that  used  for  the  Steepest  Descent  method.  That  is,  when 
an  estimate  crosses  over  the  boundary  constraint  t  it  is  then 
replaced  by  that  boundary  constraint.  The  Gauss-Newton  algorithm 
as  implemented  on  the  computer  is  presented  in  Figure  4.2. 

Scalar  Gauss-Newton  Method.  A  simplified  version  of  the 
Gauss-Newton  algorithm  was  proposed  by  Dr.  Quinn.  (QUIN83) 

The  simplification  is  recognized  when  we  let  the  Gauss-Newton 
method  update  only  one  parameter  at  a  time.  If  the  Gauss- 
Newton  algorithm  updated  only  the  parameter  qQ,  then  the 
algorithm  reduces  to 

qR(i+l)  =  q^(i)  -  (4.22) 
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Gauss-Newton  Algorithm 


with  the  starting  point  £,  evaluate  ssq 
DO  WHIIE  (counter  .LE.  max_count ) 

find  e(i,j,k)  =  (0(i, j,k)  -  y(i, j,k))/0(if j,k)  for  i=1  to  4, 

j=1  to  3, 
k=1  to  15 

FOR  n=1  to  7  DO 

with  q_test(n)  =  1.0001*q(n)  and  q_test(m)  =  q(m)  for  nj^n, 
find  e_n(n,i, j,k)=(0(i, j,k)  -  y(i, j,k))/0(i, j,k) 
dje(n,i,j,k)  =  (e_n(n,i, j,k)  -  e(i, j,k)/(1.0001*q(n) ) 
ENDFOR  n 

FOR  n  =  1  to  7  DO 
FOR  o=1  to  7  DO 
D(n)  =  0.0 
C(n,m)  =  0.0 
FOR  i=1  to  4  DO 
FOR  j=1  to  3  DO 
FOR  k=1  to  15  DO 

C(n,m)  =  C(n,m)  +  d_ e(n,i, j,k)*d_e(m,it j,k) 

D(n)  a  D(n)  +  e(it jtk)*d_e(n,i,j»k) 

ENDFOR  k 
ENDFOR  j 
ENDFOR  i 
ENDFOR  o 
ENDFOR  n 


FIGURE  4.2  Gauss-Newton  Algorithm 


-1 

evaluate  q  new  =  £  -  C  *D 
FOR  i=1  1;o  7  DO 

IF  (q_new(i)  .LT.  loverbound^ )  THEN 
q_nev(i)  =  lowerboundi 
ELSE  IF  (q__nev(i)  .IT.  upperbound^ )  THEN 
cL_nev(i)  =  upperbound^ 

ENDIF 
ENDFOR  i 

evaluate  ssq__test  for  g  new 
IF  (ssq_test  .LT.  ssq)  THEN 
£  =  £  new 
ssq  =  ssq_test 
ENDIF 

counter  +  counter  +  1 
ENDWHILE 

The  final  estimate  is  £  with  a  sum  of  s  uares  error  of  ssq 


This  algorithm  is  in  fact  equivalent  to  the  Gauss-Newton  algorithm 
when  an  elements  in  the  C  matrix  are  zero  except  for  the  diagonal 

elements. 

The  advantage  of  the  scalar  form  is  two-fold.  First,  a 
large  amount  of  computation  is  avoided  since  none  of  the  cross¬ 
terms  in  the  C  matrix  needs  to  he  computed.  In  addition,  no  matrix 
inversion  is  required.  Instead,  we  compute  the  inverse  of 
C  simply  by  inverting  each  diagonal  element.  The  Scalar  Gauss- 
Newton  algorithm  as  implemented  on  the  computer  is  presented 
in  Figure  4.3. 


Validation  of  the  Computer  Prog 


Plots  of  both  the  data  and  the  predicted  values  indicated 


absence  of  any  major  program  logic  errors  or  data  input  errors.  As 
a  check  for  further  possible  errors,  the  computer  programs'  results 
at  each  iteration  were  examined.  These  checks  indicated  that 


the  algorithms  correctly  updated  the  estimates  with 
feasible  values  that  decreased  the  value  of  the  error 


function. 


Although  the  three  optimization  algorithms  were  working 
correctly,  the  values  calculated  in  the  error  function  still 
needed  to  be  validated.  This  was  done  by  verifying  that  the 
Hexane  Model  equations  were  implemented  on  the  computer  correctly. 


Scalar  Gauss-Newton  Algorithm 

evaluate  ssq  for  the  starting  value  £(0) 
evaluate  e(i,j,k)  =  (0(i,j,k)  -  y(i, j,k))/0(i, j,k) 


for  i=1 ,4 
J-1.3 
k=1,15 


DO  WHII£  (counter  «LE.  max_count) 

FOB  n=1  to  7  DO 

q_temporary  =  q(n) 
q(n)  =  1.0001*q(n) 
evaluate  ssq__new  for  £ 

evaluate  e_new(i, j,k)  in  the  same  manner  as  was  done  for 

e(itj,k) 

q(n)  =  q_temporary 
sum_num  =  0.  ( 
sum_den  =  0. 

FOB  i  =  1  to  4  DO 
FOB  j  a  1  to  3  DO 
FOB  k  =  1  to  15  DO 

sum  num  a  sum  num  +  (e(i,j,k)  -  e  new(it j,k) )* 

e(i, j,k)/(.0001*q(n)) 

sumjden  =  sum_den  +  (e(i,j,k)  -  e__new(i,  j,k))  / 

(,0001*q(n))2 

ENDFOR  k 
ENDFOR  j 
ENDFOR  i 


FIGURE  4.3  Scalar  Gauss-Newton  Algorithm 
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q_nev(n)  =  s  um_num/ suo_de  n 
IF  (q__new(n)  .LT.  loverboundQ)  THEN 
q_new(n)  =  loverbound^ 

ELSE  IF  (q_new(n)  .GT.  upperbound^)  THEN 
q_nev(n)  =  upperbound^ 

ENDIF 
ENDFOR  n 

evaluate  ssq_test  for  q  new 
IF  (ssq^test  .LT.  ssq)  THEN 
2.  =  q  new 
ssq  =  ssq__test 
ENDIF 

counter  =  counter  +  1 
ENDWHILE 

print  out  the  results:  £  and  ssq 


FIGURE  4.3  (continued) 
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To  verify  this,  Dr.  Q^inn  and  Major  Clewell  both  independently 
solved  the  Hexane  Model  equations  on  the  computer.  When 


their  results  agreed  exactly  with  those  in  the  optimization 
♦ 

routines,  it  was  concluded  that  the  equations  in  fact  were 
programmed  into  the  optimization  routines  correctly. 

Comparison  of  the  Algorithms 

Two  types  of  comparisons  were  made  for  the  three  algorithms. 

In  the  first  comparison,  the  number  of  additions  and  subtractions, 
multiplications,  and  divisions  required  by  each  of  the  algorithms 
were  counted  for  each  iteration.  In  addition,  the  number  of  times 
the  IMSL  subroutines  DGEAR  and  LINV2F  are  called  were  counted. 

Recall  that  the  subroutine  DGEAR  solves  systems  of  stiff  differential 
equations  to  determine  the  predicted  values  of  the  Hexane  Model. 

The  subroutine  LINV2F  computes  the  inverse  of  a  matrix.  The 
list  of  the  number  of  function  evaluations  is  given  in  Table  4.1. 


TABLE  4.1  The  number  of  function  evaluations  for 
_ each  iteration  of  each  algorithm _ 


#  Additions/ 

#  Subtractions  #  Multiplications  #  Divisions  DGEAR  LINV2F 


Steepest 

Descent 

14 

21 

7 

8 

— 

Gauss-Newton 

35,672 

35,672 

17,640 

8 

1 

Scalar 

Gauss-Newton 

3,794 

3,787 

1,267 

8 

We  should  note  that  the  Steepest  Descent  computer  program 
does  not  always  pass  through  the  same  set  of  calculations  at  each 
iteration.  The  Steepest  Descent  algorithm  calculates  the  gradient 
▼actor  during  the  first  iteration.  However,  it  does  not  recompute 
the  gradient  vector  at  each  iteration.  If  the  current  gradient 
vector  fails  to  provide  updates  with  a  lower  sura  of  squares  value, 
even  when  the  step  size  J  has  been  decreased,  then  the  gradient 
vector  is  recomputed  in  the  next  iteration.  Table  4.1  lists 
the  number  of  function  evaluations  for  the  case  of  the  Steepest 
Descent  algorithm  computing  the  gradient  vector.  If  the  Steepest 
Descent  algorithm  did  not  evaluate  the  gradient  this  iteration, 
then  there  would  be  no  divisions  and  only  one  call  to  DGEAR. 

In  addition,  the  number  of  subtractions  and  the  number  of 
multiplications  each  equals  seven. 

The  comparison  in  Table  4.1  shows  that  the  Gauss-Newton 
algorithm  has  many  more  function  evaluations  than  both  the 
Steepest  Descent  algorithm  and  the  Scalar  Gauss-Newton 
algorithm.  On  a  smaller  scale,  the  Scalar  Gauss-Newton 
algorithm  has  more  function  evaluations  than  the  Steepest 
Descent  algorithm.  The  optimization  strategy  to  be  employed 
combines  the  Steepest  Descent  algorithm  and  only  one  of  the 
Gauss-Newton  algorithms.  Thus,  of  primary  concern  is  the 
comparison  between  the  two  Gauss-Newton  algorithms  since  only 
one  of  them  will  be  used.  A  second  comparison  was  made  to 


over  the  Scalar  Gauss-Newton  algorithm  is  worth  the  convergence 
speed  it  provides. 

In  the  second  comparison,  the  algorithms'  performances 
were  determined.  Each  algorithm  was  timed  for  speed 
in  executing  a  given  number  of  iterations.  Also,  the 
sum  of  squares  error  was  tabulated.  The  results  are  shown 
in  Tables  4.2  and  4.3. 


TABLE  4.2 

Time  (secs) 

Steepest 

Gauss- 

Scalar 

#  of  iterations 

Descent 

Newton 

Gauss-Newton 

1 

35-106 

193.132 

55.786 

2 

38.479 

388.954 

101.231 

#  3 

41.303 

580.838 

146.347 

5 

43.686 

967.780 

236.901 

10 

44.909 

1912.042 

459.229 

Hesults  for  the  two  Gauss-Newton  methods  were  unexpected. 
The  Scalar  Gauss-Newton  algorithm  seemed  to  achieve  slightly 
smaller  sum  of  squares  error  values  than  the  full  Gauss-Newton 
algorithm.  The  opposite  was  expected.  In  addition,  the  full 


Gauss-Newton  aethod  took  much  more  time  to  execute  each  iteration, 


This  seemingly  poorer  performance  of  the  full  Gauss-Newton 
algorithm  perhaps  can  be  attributed  to  the  numerical  errors 
that  may  be  introduced  with  the  greater  number  of  computations 
it  performs.  Based  on  the  results,  the  full  Gauss-Newton 
algorithm  was  excluded  from  any  further  optimization  runs. 


TABLE  4.3  Sum  of  squares  error 


t 


#  of  iterations 

Steepest 

Descent 

Gauss- 

Newton 

Scalar 

Gauss-Newton 

1 

32.171 

26.772 

25.361 

2 

32.171 

26.035 

24.807 

3 

32.171 

25.663 

24.513 

• 

5 

30  659 

24.970 

24.225 

10 

28.070 

24.574 

24.022 

One  further  comment  should  be  made  about  the  algorithms. 
Little  effort  was  made  to  reduce  the  calculations  involved  in 
each  algorithm  to  the  minimum.  After  the  results  shown  in 
Tables  4.2  and  4.3  were  gathered,  the  author  has  noted  that 
the  amount  of  function  evaluations  in  the  full  Gauss-Newton 
method  can  be  reduced.  However,  in  the  interest  of  time,  the 
reduced  amount  of  function  evaluations  could  not  be  implemented 


V.  Estimation  of  the  Parameters 

Once  the  optimization  routines  were  implemented  on  the 
computer,  the  estimation  process  was  initiated.  This  process 
was  iterative  in  that  estimation  results  were  frequently  passed 
to  AFAMRL  after  which  they  returned  with  further  modifications 
to  the  model.  For  parameter  estimation  of  the  Hexane  Model, 
preliminary  results  were  first  examined  by  Dr.  Anderson  and 
Major  Clewell  to  ensure  that  they  made  physical  sense.  Any 
discrepancies  that  the  researchers  noted  were  then  corrected. 
Next,  each  of  the  algorithms*  performance  was  studied  to  see 
which  step  sizes  appeared  to  be  the  most  efficient.  Third, 
several  starting  values  were  selected  as  part  of  a  strategy 
for  seeking  the  global  minimum.  Finally,  the  estimation  results 
were  examined  by  AFAMRL  before  acceptance. 


Results  for  Model  Validation 


Some  computer  runs  were  initially  made  using  starting 


values  recommended  by  Dr.  Anderson  and  Major  Clewell.  These 
routines  were  sent  to  AFAMRL  to  be  critiqued.  AFAMRL' s 


recommendations  were  then  incorporated  into  the  routines  before 
a  large  amount  of  effort  was  expended  into  estimating  the 
parameters.  As  it  turned  out,  the  data  plots  proved  invaluable 
since  AFAMRL  could  judge  for  themselves  the  performance  of  the 
optimization  routines.  The  plots  also  allowed  AFAMRL  to  quickly 
detect  some  discrepancies  between  the  model  and  the  data.  In 
light  of  these  discrepancies,  Dr.  Anderson  changed  the  following 
parameters  in  the  Hexane  Model  to 

FB  =  .8 

&  =  200. 

PR  =  11.5 
PS  =  3-62 


These  changes  provided  a  better  fit  between  the  data  and  the 
predicted  value  of  concentration  levels  of  hexane  in  the  blood. 

Step  Sizes  for  Efficient  Algorithm  Implementation 

In  the  Steepest  Descent  algorithm,  the  sum  of  squares 
error  was  examined  for  various  values  of  the  step  size  J. 

The  results  are  shown  in  Table  5.1.  It  appears  from  Table  5.1 
that  for  the  Hexane  Model,  the  Steepest  Descent  method  works 
most  efficiently  with  the  step  size  J  =  .001.  This  result, 


however,  may  not  apply  for  different  starting  values,  different 


Number  of 
Iterations 


TABLE  5.1  Sum  of  Squares  Error  for  the  Steepest 
Descent  Algorithm 


Step  Size 


.00001 

.0001 

.001 

.01 

.1 

37.306 

36.449 

34.703 

37.416 

37.416 

37.199 

35.735 

34.703 

34.703 

37.416 

36.397 

34.649 

32.525 

33.607 

34.703 

Next,  the  sum  of  squares  error  was  examined  in  the  Hexane 

Model  for  the  various  normalized  step  sizes  Aq^/q^  used  in 

the  Scalar  Gauss-Newton  method.  Recall  that  Aq  was  used  as 

n 

part  of  the  finite  difference  method  for  approximating  the 
derivatives.  The  results  are  shown  in  Table  5.2. 


TABLE  5.2  Sum  of  Squares  Error  for  the  Scalar 
Gauss-Newton  Algorithm 


Number  of 
Iterations 


Normalized  Step  Size  4qn/qQ 


The  normalized  step  size  Aqn/qQ  -  *00001  should  provide  the  best 
approximation  to  the  derivatives.  However,  Table  5«2  shows  that  the  value 
Aqn/qn  =  .001  might  be  used  without  degrading  performance.  Again,  we 
note  that  these  results  may  vary  with  the  choice  of  starting  values, 
error  function,  or  pharmacokinetic  model. 

i 

Selecting  the  Initial  Values 

After  obtaining  these  preliminary  results  about  the  Hexane  Model 
and  the  optimization  algorithms,  a  procedure  to  obtain  the 
global  optimum  was  started.  First,  an  initial  set 
of  starting  values  were  selected.  In  this  case,  the  selected 
values  were  the  upper  and  lower  bounds  of  the  unknown  parameters. 

This  gives  a  total  of  128  possible  starting  points.  The  option 
of  using  a  larger  number  of  starting  points  was  considered, 
but  it  was  noted  that  such  an  increase  was  too  large  to  handle 
in  the  time  allocated  for  this  research.  For  example,  if  we 
had  selected  3  values  instead  of  2  values  for  each  of  the  unknown 
parameters,  then  a  total  of  2,1 87  starting  points  would  need 
to  be  evaluated.  However,  it  was  concluded  that  examining 
2,187  possible  starting  points  was  too  large  a  task  for  this 
thesis  effort.  Fortunately,  AFAMR1  agreed  that  the  128  starting 
points  would  be  sufficient  to  explore  the  possibilities. 

Hence,  only  the  128  possible  starting  points  were  selected  for 


further  analysis. 

For  each  of  the  123  possible  starting  points,  the  suin  of 
squares  error  was  computed.  The  results  were  then  scanned 
so  that  points  with  a  sum  of  squares  error  below  45  could  be 
selected  as  starting  points.  The  other  starting  points  were 
discard.  The  12  points  with  a  sum  of  squares  error  below  45  are 
shown  in  Table  5*3* 


TABLE  5«3  Selected  Starting  Values 


t 

KFO  VMAX2  KM2  VMAX3  KM3  KM2I 


In  addition  to  the  12  starting  points  presented  in  Table  5*3 
AFAMRL  provided  starting  values  based  upon  what  they  felt 
would  be  the  best  estimate  of  the  parameters*  Their  starting 
values  are  shown  in  Table  5*^*  Altogether)  13  starting  points 
were  used  by  the  optimization  algorithms* 


TABLE  5*4  Starting  Values  Recommended  by 
AFAMRL 


Each  of  the  12  points  in  Table  5*3  were  first  used  as 
starting  values  for  the  Steepest  Descent  algorithm.  The  Steepest 
Descent  algorithm  then  updates  the  initial  values  3  times. 

These  few  iterations  are  intended  to  bring  the  initial  values 
closer  to  a  local  minimum.  The  resulting  estimates  are  then 
used  for  further  optimization  with  the  Scalar  Gauss-Newton  method 
Since  the  initial  values  are  relatively  close  to  a  local  minimum 
after  passing  through  the  Steepest  Descent  algorithmT  the 


1000.000  .315 


In  addition  to  the  12  starting  points  of  Table  5.3,  the 
starting  values  provided  by  AFAMRL  underwent  3  iterations  of  the 
Steepest  Descent  algorithm  and  then  15  iterations  of  the  Scalar 
Gauss-Newton  algorithm.  The  results  for  this  case  are  shown  in 
Table  5.6. 


TABUS  5.6  Results  from  the  starting  values  provided 
provided  by  AFAMRL  in  Table  5.4 


KFO  WMAX2  KM2 
2.4  .4025  .5317 


VMAX3  KM3  KM2I 

2.8177  2.5132  1000.0 


K 


2931 


SSQ 

28.428 


Of  the  13  results  presented  in  Tables  5.5  and  5.6,  the 
starting  point  provided  by  AFAMRL  showed  the  lowest  value 
in  its  sum' of  squares  error.  The  results  shown  in  Table  5.6 
were  thus  forwarded  to  AFAMRL  for  approval.  The  major  factor 
influencing  AFAMRL1 s  decision  to  accept  the  estimate  was  not 
necessarily  the  low  value  for  the  sum  of  squares  error. 

Rather,  the  plots  provided  a  clear  picture  of  the  performance 
of  the  estimates  to  AFAMRL.  The  estimates  were  thus  accepted 
since  the  plots  showed  that  for  the  Hexane  Model,  the  predicted 
values  match  the  data  closely  in  rainy  of  the  plots.  The  plots 
are  shown  in  Figure  5.1. 


A  final  comment  should  be  made  about  the  results.  With  the 
many  updates  that  AFAMRL  made  to  the  model,  the  final  error  function 
in  the  algorithms  did  not  exactly  represent  the  Hexane  Model  used 
in  this  thesis.  However,  as  it  turns  out,  the  error  function  still 
measured  the  agreement  between  the  predicted  values  and  the 
experimental  data.  The  estimates,  although  perhaps  not  global, 
served  their  purpose  by  indicating  to  the  scientists  the 
underlying  processes  involved  with  the  hexane  in  the  mammal. 

The  results  met  AFAMRL* s  approval. 


Concentration  of  hexane  in  the  blood:  Run  ft' 


FIGURE  5.1(d)  Concentration  of  raethyl-n-butylketone  in  body  water:  Run  #2 


Concentration  of  2,5-bexanedione  in  body  water 
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(h)  Concentration  of  2, 5-hexane di one  in  body  water:  Run  #5 


FIGURE  5*1(1)  Concentration 


Concentration 

(ag/L) 


FIGURE  5*1 (j)  Concentration  of  methyl -n-butylketone  in  body  water:  Run  #4 


Concentration  of  2,5-hexanedione  in  body  water:  Run  #4 


VI.  Sensitivity  Analysis 


I 

For  the  estiaates  of  the  unknown  parameters  obtained  in 
Chapter  V,  a  sensitivity  analysis  was  conducted.  The  analysis 
was  intended  to  aid  the  researchers  at  AFAMRL  in  determining 
which  set  of  the  estimated  parameters  had  a  large  impact  on 
the  error  function.  One  of  the  purposes  for  determining  this 
is  to  help  AFAMRL  in  designing  future  experiments.  By  understanding 
which  parameters  affected  the  error  function  heavily,  the 
researchers  could  then  experiment  in  an  attempt  to  estimate 
precisely  the  more  influential  parameters.  In  the  meantime, 
the  least  influential  parameters  may  be  eliminated  from  further 
estimation. 

In  this  chapter,  a  post-experimental  sensitivity  analysis 
was  conducted  first.  This  type  of  analysis  is  typically  conducted 
after  the  experimentation  and  parameter  estimation.  From  it 
come  clues  about  how  a  change  in  the  estimates  would  affect 
the  error  function.  In  essence,  post-experimental  sensitivity 
analysis  indicates  the  degree  to  which  each  parameter  is 
estimated  precisely. 


The  second  type  of  analysis  involves  the  examination  of 
sensitivity  functions.  This  type  of  analysis  may  be  used  before 
experimentation*  For  this  thesis t  the  analysis  of  sensitivity 
functions  was  relatively  brief.  The  intent  of  this  analysis  was 
to  introduce  the  researchers  to  the  possibilities  of  experimental 
design  of  dynamical  systems*  Although  brief*  the  discussion 
of  this  type  of  analysis  should  provide  the  researchers  a 
useful  starting  point  to  designing  experiments  with  sensitivity 
functions. 

Post-Experimental  Sensitivity  Analysis 

For  post-experimental  sensitivity  analysis*  the  approach 
used  is  similar  to  that  recommended  by  Lehman  and  Stark. 

To  begin  with*  each  of  the  values  were  evaluated  at  a  factor  of 
-  .00001  away  from  their  estimates.  The  sum  of  squares  error 
for  each  of  these  perturbed  values  is  shown  in  Table  6.1.  Next* 
the  factor  was  increased  by  a  multiple  of  3*  Thus*  a  second  set 
of  perturbed  values  had  their  sum  of  squares  error  computed  at 
a  factor  of  -  .00003  away  from  their  estimates.  This  process 
of  increasing  the  perturbation  factor  by  a  multiple  of  3  was 
continued  several  times.  For  each  perturbed  value*  the  sum  of 
squares  error  was  computed.  The  results  are  displayed  in 
Table  6.1. 
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TAME  6.1(a)  Sensitivity  Analysis  for  KFO 


Factor  A 

KF0*(1  -  A) 

SSQ 

KF0*(1  +4) 

ssq 

.00001 

2.399976 

28.430 

2.400024 

28.- 

.00003 

2.399928 

28.430 

2.400072 

28.j 

.00009 

2.399784 

28.430 

2.400216 

28J 

.00027 

2.399352 

28.430 

2.400648 

28.' 

.00081 

2.398056 

28.430 

2.401944 

28. ; 

.00243 

2.394168 

28.429 

2.405832 

28j 

.00729 

2.382504 

28.427 

2.417496 

28. J 

.02187 

2.347512 

28.411 

2.452488 

28J 

.06561 

2.242536 

28.411 

2.557464 

28J 

.19683 

1.927608 

28.455 

2.872392 

28.1 

.59049 

.982824 

29.092 

3.817176 

30. 

TAME  6.1(b)  Sensitivity  Analysis  for  VMAX2 


Factor  4 

VMAX2*(1  -  A) 

SSQ 

VMAX2*(1  +a) 

SSQ 

.00001 

.4025256 

28.430 

.4025336 

28. 

.00003 

.4025175 

28.430 

.4025417 

28. 

.00009 

.4024934 

28.430 

.4025658 

28. 

.00027 

.4024209 

28.430 

.4026383 

28. 

.00081 

.4022036 

28.430 

.4028556 

28. 

.00243 

.4015515 

28.429 

.4035078 

28. 

.00729 

.3995952 

28.427 

.4054641 

28. 

.02187 

.3937263 

28.425 

.4113329 

28. 

.06561 

.3761196 

28.451 

.4289396 

28.1 

.19683 

.3232998 

28.835 

.4817595 

28.' 

.59049 

.1648399 

33.402 

.6402193 

32.- 
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TABLE  6.1(c)  Sensitivity  Analysis  for  KM2 


Factor  A 

.00001 

.00003 

.00009 

.00027 

.00081 

.00243 

.00729 

.02187 

.06561 

.19683 

.59049 


KM2*(1  -  A) 

SSQ 

KM2*(1+A) 

.5317015 

28.430 

.5317121 

.5316908 

28.430 

.5317227 

.5316589 

28.430 

.5317546 

.5315632 

28.430 

.5318503 

.5312761 

28.430 

.5321375 

.5304147 

28.430 

.5329988 

.5278306 

28.431 

.5355829 

.5200783 

28.433 

.5433352 

.4968215 

28.451 

.5665921 

.4270509 

28.652 

.6363626 

.2177392 

32.441 

.8456743 

28.430 

28.430 

28.430 

28.430 

28.430 

28.430 

28.430 

28.432 

28.447 

28.564 

29.266 


TABLE  6.1(d)  Sensitivity  Analysis  for  VMAX3 


Factor  A 

.00001 

.00003 

.00009 

.00027 

.00081 

.00243 

.00729 

.02187 

.06561 

.19683 

.59049 


VMAX3*(1  - 

2.817675 
2.817619 
2.816943 
2.816943 
2.815422 
2.810857 
2.797163 
2.756081 
2.632835 
2.26 3095 
1.153878 


)  SSQ 

VMAX3*(1  +A) 

SSQ 

28.430 

2.817732 

28.430 

28.430 

2.817789 

28.430 

28.430 

2.818465 

28.431 

28.430 

2.818465 

28.430 

28.430 

2.819986 

28.431 

28.429 

2.824551 

28.433 

28.432 

2.838245 

28.443 

28.479 

2.879327 

28.507 

29.047 

3.002574 

28.957 

38.166 

3.372313 

31.509 

1199.806 

4.481530 

39.159 

.V*."  ■«  \ 


TA3LE  6.1(e)  Sensitivity  Analysis  for  KM3 


Factor  A 

KM3*(1  -4) 

SSQ 

KM3*(1  +4) 

SSQ 

.00001 

2.513165 

28.430 

2.513216 

28.430 

.00003 

2.513115 

28.430 

2.513266 

28.430 

.00009 

2.512964 

28.430 

2.513417 

28.430 

.00027 

2.512512 

28.430 

2.513869 

28.430 

.00081 

2.511155 

28.430 

2.515226 

28.430 

.00243 

2.507083 

28.430 

2.519297 

28.430 

.00739 

2.494869 

28.430 

2.531512 

28.431 

.02187 

2.458226 

28.436 

2.568154 

28.439 

.06561 

2.343300 

28.495 

2.678081 

28.500 

.19683 

2.018519 

29.095 

3.007862 

28.992 

.59049 

1.029177 

36.640 

3.997204 

32.695 

TABLE  6.1(f)  Sensitivity  Analysis  for  KM2I 


Factor  A 

KM2I*(1  -4) 

SSQ 

KM2I*(1  +4) 

SSQ 

.00001 

999.99 

28.430 

1000.01 

28.430 

.00003 

999.97 

28.430 

1000.03 

28.430 

.00009 

999.91 

28.430 

1000.09 

28.430 

.00027 

999.73 

28.430 

1000.27 

28.430 

.00081 

999.19 

28.430 

1000.81 

28.430 

.00243 

997.57 

28.430 

1002.43 

28.430 

.00729 

992.71 

28.430 

1007.29 

28.430 

.02187 

978.13 

29.430 

1021.87 

28.430 

.06561 

934.39 

28.430 

1065.61 

28.430 

.19683 

803.17 

28.431 

1196.83 

28.430 

.59049 

><09.5 1 

28.436 

1590.49 

28.429 

TABLE  6.1(g)  Sensitivity  Analysis  for  K 


Factor  A 

K*(1  -  A) 

.00001 

.2930817 

.00003 

.2930758 

.00009 

.2930582 

.00027 

.2930055 

.00081 

.2928472 

.00243 

.2923724 

.00729 

.2909480 

.02187 

.2866748 

.06561 

.2738553 

.19683 

.2353968 

.59049 

.1200211 

SSQ 

K*(1  +A) 

28.430 

.2930875 

28.430 

.2930934 

28.430 

.2931110 

28.430 

.2931637 

28.430 

.2933220 

28.430 

.2937970 

28.432 

.2952212 

28.444 

.2994944 

28.560 

.3123139 

30.061 

.3507724 

84.209 

.4661481 

From  Table  6.1,  it  appears  that,  in  general,  very  small  changes 


in  all  of  the  parameters  do  not  affect  the  error  function.  The 
parameter  VMAX3,  however,  affects  the  error  function  much  more 
than  all  other  parameters  at  most  of  the  perturbation  values. 

At  the  perturbation  factor  of  .59049,  VMAX3  has  the  largest 
sum  of  squares  error  of  1199«8o6.  Comparing  the  other  error 
values  in  VMAX3  with  that  of  the  other  parameters,  we  conclude 
that  the  estimate  for  VMAX3  must  be  known  as  closely  as  possible 
in  order  that  the  model  closely  fit  the  data. 

The  parameter  KM2I,  on  the  other  hand,  does  not  appear  to 
affect  the  error  function  at  all.  This  implies  that  the 
estimate  for  KM2I  may  not  need  to  be  precise  at  all. 

AFAMRL  may  consider  leaving  KM2I  as  a  known  constant  in  future 
experiments  since  its  value  may  have  little  effect  on  the 
final  estimates. 


Sensitivity  Functions  for  Experimental  Design 

The  use  of  sensitivity  functions  for  experimental  design 
may  be  valuable  for  AFAMRL.  The  researchers  are  concerned 
about  how  they  might  conduct  future  experiments  with  the  Hexane 
Model.  In  pharmacokinetic  applications,  it  appears  that  the 
number  of  data  points  from  an  experiment  is  not  necessarily  large. 


In  addition,  the  length  of  the  experiment  may  be  constrained  by 
time  limits.  With  such  limitations  on  the  experiment,  a  careful 
plan  for  experimentation  becomes  important.  Hence,  the  analysis 
of  sensitivity  functions  before  experimentation  is  appropriate 
for  .pharmacokinetic  applications. 

The  analysis  of  sensitivity  functions  covered  in  this  section 

is  relatively  brief.  The  topic  itself  can  be  involved  enough 

to  warrant  a  thesis  effort  solely  devoted  to  experimental  design 

topics  for  dynamical  systems.  The  intent  of  this  section  is  to 

introduce  AFAMRL  and  other  readers  to  the  applications  of  sensitivity 

functions.  To  do  this,  we  first  consider  the  problem  of  weighting 

the  sensitivity  functions.  This  problem  has  not  been  resolved 

in  this  thesis  effort,  but  it  can  be  the  starting  point  of  another 

thesis.  Next,  some  sensitivity  functions  are  considered  for  the 

experiment  which  has  already  been  conducted  on  the  Hexane  Model. 

By  doing  this,  we  can  compare  what  we  already  know  from  the 

post-experimental  sensitivity  analysis  with  the  results  of  the 

sensitivity  functions.  Next,  a  small  example  of  experimental 

design  is  presented.  Finally,  a  brief  comment  is  made  about 

« 

optimal  input  design. 

The  Proper  Weighting  of  Sensitivity  Functions.  Chapter  III 


presented  some  concepts  about  using  sensitivity  functions  for 
experimental  design.  However,  it  is  still  not  clear  from  the 


literature  exactly  how  we  might  go  about  using  these  sensitivity 
functions*  For  example,  Lehman  and  Stark  suggest  plotting  the 
sensitivity  functions 

in  order  to  get  a  feeling  for  which  parameters  have  the  most 
influence  on  each  of  the  measurements  y^*  Lehman  and  Stark 
then  recommend  that  the  sensitivity  functions  be  normalized  by 
the  nominal  parameter  q&.  This  way,  the  plots  of  the  magnitudes 
of  the  sensitivity  functions  can  be  directly  compared  with  each 
other.  Then  the  plots  with  the  larger  magnitudes  identify  the 
parameters  that  are  likely  to  have  a  large  impact  on  the  error 
function*  The  plots  also  show  the  points  in  time  when  cm  experimental 
observation  is  likely  to  have  a  large  effect  in  estimating  the 
parameter  under  consideration*  (LEHM81 :  114-117) 

As  reported  in  Chapter  HI,  Beck  and  Arnold  recommend  that 
for  single  response  'cases,  the  sensitivity  functions  be  normalized 
in  the  same  manner  as  that  recommended  by  Lehman  and  Stark.  However, 
for  the  multiple  response  case,  Beck  and  Arnold  add  that  the  sensitivity 
function  be  further  normalized  by  the  standard  deviation  of  the  measuring 
instruments.  From  Beck  and  Arnold's  discussion,  it  appears  that 
in  order  to  compare  the  relative  magnitudes  of  the  sensitivity 


functions,  we  should  normalize  the  functions  with  both  the 


nominal  values  q^  and  the  standard  deviations  of  the  measurements 
<^.  Then,  if  the  standard  deviations  happen  to  be  equal  for 
all  measurements,  we  only  need  to  normalize  with  the  nominal 
parameters  as  suggested  by  Lehman  and  Stark* 

Although  Arnold  and  Beck's  suggestions  on  normalizing  the 
sensitivity  functions  appear  to  be  clear,  there  is  confusion 
in  its  applicability.  It  seems  that  Beck  and  Arnold  were  making 
their  suggestions  of  normalizing  the  sensitivities  for  the  case 
of  the  unweighted  least  squares  error  function*  (BECK77:  346-349) 
No  comment  is  made  concerning  the  impact  of  estimating  with  the 
weighted  least  squares  error  function.  Lehman  and  Stark  do 
not  specify  at  all  what  their  error  criterion  is.  It  may  be 
possible,  however,  that  they  too  were  announcing  techniques  that 
assumed  an  unweighted  least  squares  criterion.  At  any  rate, 
there  are  uncertainties  in  the  author's  mind  as  to  how  we 
might  properly  weight  the  sensitivity  functions  in  order  to 
'  directly  compare  their  relative  magnitudes.  It  is  not  well  known 
by  the  author  how  the  choice  of  the  error  function  relates  to 
the  proper  weights. 

For  the  case  of  the  Open  Chamber  Hexane  Model,  the  standard 
deviations  of  the  measurements  are  not  known.  Thus,  we  could 
not  properly  normalize  the  sensitivities  by  their  standard 
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deviations  even  if  we  wanted  to.  In  addition,  since  the  weighted 
least  squares  error  function  was  used,  we  are  not  sure  how  the 
sensitivity  functions  should  be  normalized  in  order  to  compare 
exactly  the  effects  of  the  parameters  on  the  error  function. 

Thus,  only  the  unweighted  sensitivity  functions  were  considered 
and  plotted  in  this  section.  These  uncertainties  about  weighting 
the  sensitivity  functions  most  likely  can  be  resolved  with  some 
effort.  The  effort  required,  however,  is  beyond  the  scope  of 
this  thesis  project.  Still,  the  unweighted  analysis  can  serve  us 
as  an  introduction  to  the  reader  on  how  an  experiment  can  be 
designed. 

Some  Sensitivity  Functions  of  the  Current  Experiment  on  the 
Hexane  Model.  Before  describing  the  sensitivity  plots,  we  first 
recall  that  in  the  Hexane  Model  there  were  7  unknown  parameters. 

This  fact  implies  that  7  sensitivity  functions  exist  for  every 
differential  equation  in  the  Hexane  Model.  But,  since  there  are 
only  three  measurement  models  in  the  Hexane  Model,  only  the 
sensitivity  equations  for  these  three  measurements  need  examination. 
With  three  equations,  the  total  number  of  sensitivity  functions  becomes 
21.  However,  the  21  sensitivity  equations  represent  one  experimental 
run.  Since  there  were  4  experimental  runs  each  at  different 
chamber  concentration  levels,  then  a  total  of  84  plots  of  the 


sensitivity  functions  is  needed  to  study  experimental  influences  on  the 


Model*  To  get  an  indication  of  how  such  a  sensitivity  analysis 
proceeds ,  3  plots  of  the  sensitivity  functions  are  discussed* 

i 

These  3  plots  should  be  sufficient  to  introduce  the  basic 
concepts. 

In  Figure  6.1,  three  sensitivity  functions  involving  VMAX3 
and  KM2I  are  plotted.  The  plots  were  chosen  based  upon  the  results 
of  the  post-experimental  sensitivity  analysis.  It  appeared 
that  from  the  post-experimental  sensitivity  analysis  the  parameter 
VMAX3  had  the  most  influence  on  the  error  function  whereas  the 
parameter  KM2I  had  the  least  influence.  A  further  examination  of 
the  sensitivity  functions  invloving  the  parameters  VMAX3  and  KM2I 
is  likely  to  add  insight  into  the  results  already  obtained  in 
the  post-experimental  sensitivity  analysis.  Thus,  the  choice 
was  made  to  examine  some  sensitivities  involving  these  two 
parameters. 

Figure  6.1  presents  sensitivity  functions  plotted  for 
2 , 5-hexanedione  in  body  water  for  the  first  run  of  both  VMAX3 

e 

and  KM2I.  We  see  that  VMAX3  has  a  much  larger  magnitude  than 
KM2I.  This  seems  to  support  the  results  from  the  post-experimental 
sensitivity  analysis  that  VMAX3  has  a  larger  effect  on  the  measurements 
than  KM2I.  Still,  we  cannot  say  for  sure  that  the  relative 
magnitudes  in  Figure  6.1  are  directly  comparable  unless  we  had 


weighted  them  properly. 

Another  point  of  interest  is  that  both  VMAX3  and  KM2I  each 
give  the  largest  sensitivity  values  at  the  time  of  6  hours. 

/ 

This  fact  implies  that  the  experimental  observations  collected 
from  the  concentration  levels  of  2,5-hexanedione  in  body 
water  at  the  sixth  hour  should  have  the  largest  impact  on 

a 

measuring  the  parameters  KM2I  and  VMAX3  precisely.  Again, 
without  the  proper  weightings,  we  are  not  sure  how  much  will 
be  the  relative  effects. 

A  second  comparison  is  next  made  only  for  the  parameter 
VMAX3,  but  at  different  concentration  levels.  The  plot  shows 
the  sensitivities  of  2,5-hexanedione  for  runs  #1  and  #2.  The 
sensitivity  for  run  #2  is  larger  in  magnitude  than  that  for 
run  #1.  These  two  curves  represent  the  same  sensitivity  functions, 
although  the  concentration  levels  of  the  hexane  in  the  test 
chamber  is  different.  Hence,  these  two  curves  are  comparable. 

These  plots  imply  that  VMAX3  affects  the  concentration 
level  of  2,5-hexanedione  in  body  water  in  run  #2  more  *-b»n 
the  concentration  level  in  run  it 1. 

Example  of  Using  Sensitivity  Functions  to  Design  an  Experiment 
The  four  general  topics  under  experimental  design  of  dynamical  systems 
were  stated  by  Goodwin  and  Payne  and  quoted  in  Chapter  III.  This 
section  will  briefly  describe  how  one  of  these  topics  might  be 


employed  on  the  Hexane  Model.  Specifically ,  the  topic  is  input 

design.  Attention  was  placed  on  examin 

input  to  increase  the  precision  of  the  estimates. 

In  this  section,  only  one  input  chamber  concentration  level  of 
hexane  was  used.  Furthermore,  only  two  sensitivity  functions  are 
considered.  However,  the  procedure  for  examining  the  sensitivity 
plots  in  this  case  also  apply  to  other  input  concentration  levels 
and  other  sensitivity  functions. 

The  chosen  input  concentration  level  of  hexane  in  this  section 
is  shown  in  Figure  6.2.  The  input  chosen  is  a  binary  signal  that 
alternates  between  the  two  values  0  mg/1  and  3*5  mg/1.  The 
question  we  want  to  answer  is  how  can  we  experiment  in  order  to 
estimate  the  unknown  parameters  more  precisely?  The  two 
parameters  that  we  shall  arbitrarily  consider  here  are  KM2  and 
KFO. 

Figure  6.3  shows  the  sensitivity  of  methyl-n-butylketone  in 
the  body  water  with  respect  to  KM2.  The  plot  is  basically  periodic 
as  might  be  expected  from  a  periodic  binary  input  function. 

The  highest  values  are  at  times  1  hour,  4  hours,  7  hours,  ...  . 

This  indicates  that  our  estimates  of  KM2  will  probably  be  the 
most  precise  when  the  data  is  collected  at  those  times. 

The  plot  of  the  sensitivity  of  2,3-hexanedione  in  body  water 
with  respect  to  KFO  is  shown  in  Figure  6.4.  This  plot  shows  that 


Input  chamber  concentration  level 


FIGURE  6.4  Sensitivity  function  of  2,5-hexanedione  in  body  ivater  with  respect 


in  order  to  get  relatively  high  precision  in  estimating  KFO,  we 
should  obtain  experimental  observations  of  2,5-hexanedione  in  the 
body  water  at  the  points  in  time  where  the  sensitivity  is  large. 

Some  of  the  time  points  are  7. 5  hours,  10.5  hours,  13«5  hours, 
and  16.5  hours.  We  might  also  sample  at  the  points  1.5  hours  and 
4.5  hours,  although  these  points  do  not  have  sensitivity  values 
as  large  as  the  other  points. 

Now  both  plots  involving  the  parameters  KM2  and  KFO  are 
simultaneously  compared  to  further  determine  how  we  should  conduct 
the  experiment.  From  the  plots  it  appears  that  the  most  important 
experimental  observation  point  for  estimating  KM2  is  at  the  time 
of  1  hour.  This  is  the  point  of  the  maximum  difference  between 
the  sensitivities  involving  KM2  and  KFO  with  KM2  on  the  high  side. 

At  this  point  we  should  expect  that  an  experimental  observation 
should  affect  the  precision  of  KM2  with  little  influence  from 
the  effects  of  KFO. 

By  the  same  token,  we  note  that  at  the  local  maximum  points 
of  the  sensitivity  function  involving  KFO,  the  sensitivity  function 
involving  KM2  approaches  the  value  of  zero.  At  the  points  1.5  hours 
4.5  hours,  7.5  hours,  and  so  on,  the  parameter  KFO  is  likely  to 
influence  the  experimental  observations  much  more  than  the  parameter 
KM2.  Thus,  by  gathering  data  at  these  time  points,  we  enable  our¬ 
selves  to  more  closely  resolve  the  true  value  of  KFO  during 


experimentation. 

The  brief  examination  we  made  of  two  sensitivity  plots 
shows  how  one  might  be  able  to  anticipate  to  some  degree  the 
precision  of  his  estimates.  This  examination  in  fact  can  be 
applied  to  the  many  types  of  inpat  concentration  levels  of 
interest.  Then,  from  the  examination,  the  researcher  may  be  able 
to  decide  which  input  level  best  suits  his  objectives. 

Sensitivity  Functions  for  Designing  Optimal  Inputs.  An 


extension  of  the  previously  stated  methods  of  analyzing  sensitivity 
functions  is  the  procedure  involving  input  optimization.  In 
studying  the  plots  of  the  sensitivity  functions,  an  infinite 
number  of  inputs  can  be  tried.  However,  using  calculus  of 
variations  techniques,  we  can  find  input  trajectories  that 
afford  the  most  precision  possible  with  respect  to  a 
defined  error  function.  Using  this  technique,  we  can  avoid  searching 
many  input  functions  in  order  to  select  the  one  we  desire  for 
the  experiment.  (KALA82:  225-378) 
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VII.  Conclusions  and  Recommendations 

In 'sum,  the  major  part  of  AFAMRL* 8  problem  was  successfully 
solved.  That  is,  the  unknown  parameters  of  the  Open  Chamber 
Hexane  Model  were  assigned  values  that  met  approval  of 
the  researchers  at  AFAMRL.  We  note  here  that  the  primary 
factor  in  AFAMRL’ s  acceptance  of  the  final  estimates  in  Chapter  V 
was  the  information  provided  by  the  data  plots.  It  is  recommended 
that  future  parameter  estimation  algorithms  include  computer  codes 
that  plot  the  results  when  needed. 

One  interesting  result  is  the  success  of  the  Scalar  Gauss-Newton 
algorithm  presented  in  Chapter  IV.  This  algorithm  is  simpler  to 
implement  and  faster  to  execute  than  the  Gauss-Newton  algorithm. 

In  addition  to  these  benefits,  the  Scalar  Gauss-Newton  algorithm 
also  converges  more  in  the  sum  of  squares  each  iteration  than  the 
Gauss-Newton  algorithm.  As  noted,  however,  is  the  fact  that 
little  effort  was  made  to  optimize  the  codes  so  that  they 
would  require  the  smallest  amount  of  time  possible  to  perform 
each  iteration.  It  is  recommend  that  these  three  algorithms 
in  addition  to  any  other  algorithms  still  be  considered  in 


future  applications  of  pharmacokinetic  models  until  they  have 
been  fully  evaluated  in  terms  of  their  performance. 

Another  point  of  interest  brought  out  in  this  thesis  is 
the  use  of  sensitivity  functions.  In  relation  to  the  literature 
that  exists,  the  discussion  on  sensitivity  functions  was  brief. 
However,  the  discussion  was  extensive  enough  to  indicate  to 
AFAMRL  how  they  may  resolve  some  of  their  problems  in  designing 
experiments  that  afford  more  precision  for  their  estimates. 

The  author  feels  that  further  research  and  study  is  needed  on 
this  subject  before  a  high  degree  of  confidence  can  be  achieved 
in  using  sensitivity  functions.  However,  the  power  offered 
by  sensitivity  functions  in  designing  experiments  is  sufficiently 
to  warrant  such  a  study. 

To  the  author's  knowledge,  the  only  place  where  the  equations 
of  the  Open  Chamber  Hexane  Model  are  published  is  in  this  thesis. 
Originally,  the  equations  were  presented  in  a  much  different  form. 

A  great  deal  of  effort  was  made  towards  putting  the  model  equations 
into  the  form  of  a  system  of  first  order  differential  equations. 

It  is  recommended  that  the  researchers  put  future  model  equations 
into  the  from  of  differential  equations.  This  will  help  to  avoid 
the  errors  that  occur  when  someone  else  tries  to  formulate 
the  equations  in  order  to  implement  them  on  the  computer.  In 
addition,  much  time  will  be  saved  if  the  equations  given  by  AFAMRL 


are  in  the  form  that  can  be  implemented  on  the  computer  immediately 
Finally ,  the  author  gained  most  of  his  conficence  in  the 
final  estimates  after  examining  many  values  closely  surrounding 
the  estimates*  It  is  recommended  that  AFAMRL  devote  some  time 
towards  exploring  the  values  around  the  estimates  that  they 
obtain*  The  estimates  are  not  guaranteed  to  be  global,  and 
some  type  of  scheme  should  be  employed  to  indicate  how  close 
the  estimates  might  be  to  the  global.  Only  after  examining 
the  many  values  for  the  estimates  can  we  be  sure  that  our  estimates 
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